n 


tfwW 


AD-A230  578 


in  lb  i  III 


TECHNICAL  REPORT  EL-90-11 


(£> 


MS  Army  Corps 
of  Engineers 


THREE-DIMENSIONAL,  LAGRANGIAN  RESIDUAL 
TRANSPORT  COMPUTED  FROM  AN  INTRATIDAL 
HYDRODYNAMIC  MODEL 


Mark  S.  Dortch 

Environmental  Laboratory 

DEPARTMENT  OF  THE  ARMY 
Waterways  Experiment  Station,  Corps  of  Engineers 
3909  Halls  Ferry  Road,  Vicksburg,  Mississippi  39180-6199 


DTI 


V-M:*CTE  j|% 
JAN  08 1991  1  | 

|lj§?  ILJi 


V$7  /O' 

November  1990 
Final  Report 


Approved  for  Public  Release;  Distribution  Unlimited 


Prepared  for  US  Army  Engineer  District,  Baltimore 
Baltimore ,  Maryland  21 203- 1715 


G> 


_ Unclassified _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 


la.  REPORT  SECURITY  CLASSIFICATION 
Unclassified 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 
Technical  Report  EL-90- 11 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

-  USAEWES  ' 


■  i  V*  •  «•)  il  ifsl  IH  IH*  •la)  mf. 


6c.  ADDRESS  (City,  State,  and  ZIP  Code) 
3909  Halls  Ferry  Road 
.Vicksburg,  MS  39180-6199 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 

USAED,  Baltimore 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 
Baltimore,  MD  21203-1715 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


3.  DISTRIBUTION /AVAILABILITY  OF  REPORT 
Approved  for  public  release;  distribution 
unlimited. 


S.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


7a.  NAME  OF  MONITORING  ORGANIZATION 


7b.  ADDRESS  (City,  State,  and  ZIP  Code) 


8b.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(if  applicable) 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROJECT 

TASK 

NO. 

NO. 

II.  TITLE  (Include  Security  Classification) 

Three-Dimensional,  Lagrangian  Residual  Transport  Computed  from  an  Intratidal 
Hydrodynamic  Model 


12.  PERSONAL  AUTHOR(S) 
Dortch,  Mark  S. 


13a.  TYPE  OF  REPORT  13b.  TIME  COVERED  14.  DATE  OF  REPORT  (Year,  Month,  Day)  15.  PAGE  COUNT 

Final  report  from _ TO _  November  1990  280 


16.  SUPPLEMENTARY  NOTATION 

Available  from  National  Technical  Information  Service,  5285  Port  Royal  Road,  Springfield, 
VA  22161. 


17.  COSAT1  CODES  |  18.  SUBJECT  TERMS  (Continue  on  reverie  if  necessary  and  identify  by  block  number) 

FIELD  |  GROUP  |  SU8-GROUP  | 

See  reverse. 


19*  ABSTRACT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number), ' 

_ j*.  A  procedure  was  developed  for  computing  three-dimensional.  Lagrangian  residual  cir¬ 
culation  from  an  intratidal^i.e.  contains  tidal  fluctuations^hydrodynamic  model  for 
driving  an  intertidal^!,  e.  tidal  fluctuations  are  rsmoved^transport  model.  A  three- 
dimensional,  finite  difference,  hydrodynamic  model,  CH3D,  which  uses  boundary-fitted  coor¬ 
dinates  in  planform  and  vertical  cartesian  coordinates,  was  indirectly  coupled  to  a  water 
quality  transport  model  that  uses  an  integrated  compartment  solution.  The  coupling>of  the 


two  models5was  accomplished  through  development  of  an  interface  processor  implemented 
within  the  hydrodynamic  model.  The  processor  converts  nondimensional,  contravariant  veloc¬ 
ities  in  transformed  coordinates  to  dimensional,  physical  flows  for  the  transport  model. 

The  sum  of  Eulerian  residual  velocities  and  Stokes'  drift  was  used  as  a  first-order 
approximation  for  the  Lagrangian  residual  .currents.  •Xh^Stokes'  drift  approximates 

^(Continued) 


720.  DISTRIBUTION/AVAILABILITY  OF  ABSTRACT 
[3  UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT 
22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 


DD  Form  1473,  JUN  86 


□  DTIC  USERS 


21.  ABSTRACT  SECURITY  CLASSIFICATION 


22b.  TELEPHONE  (Include. Area  Code)  22c.  OFFICE  SYMBOL 


Previous  editions  are  obsolete.  SECURITY  CLASSIFICATION  OF  THIS  PAGE 

Unclassified 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 


18.  SUBJECT  TERMS  (Continued) . 

Chesapeake  Bay 
CH3D 

Eulerian  residuals 
Hydrodynamic  and  water  quality 
model  interfacing 

Intratidal  and  intertidal  transport 


Lagrangian  residual  circulation 

Salinity 

Stokes'  drift 

Three-dimensional  numerical  model 
Tide-induced  currents 


19.  ABSTRACT  (Continued). 

•residual  currents  induced  by  •tl^f'nonlinear  interactions  of  tidal  currents  and  repre¬ 
sents  the  net  drift  experienced  by  a  particle  passing  through  a  spatially  varying  veloc¬ 
ity  field  in  an  oscillating  flow.  A  Stokes'  drift  formulation  that  guarantees  mass 
conservation  was  implemented  within  the  Interface  processor  so  that  intertidalj(i.e. 
tidally  averagedj^ydrodynamic  Information  could  be  processed~and  output  as  tire  intra¬ 
tidal  hydrodynamic  model  is  executing.  This  information  is-/s"ubsequently&used  to  drive 
intertidal  mass  transport.  K-Ciju’tn  <■)& ; 

Computed  residual  velocities  were  compared  against  a  two-dimensional  (vertical- 
longitudinal)  analytical  result  to  ensure  proper  implementation.  Other  than  some  adjust¬ 
ment  of  eddy  viscosity  to  account  for  the  effects  of  numerical  dampening,  the  computed 
residual  currents  compared  favorably  with  the  analytical  result. 

The  methodology  was  evaluated  through  an  application  to  Chesapeake  Bay.  Salinity 
computed  by  the  two  models  and  observed  data  were  compared  for  the  entire  year  1985. 
Salinity  computed  with  the  intertidal  transport  model  (i.e.  with  Lagrangian  residuals) 
showed  good  agreement  with  observed  salinity  and  that  computed  with  intratidal  transport. 

The  developed  procedures  provide  a  practical  means  of  making  long-term,  three- 
dimensional  transport  computations  for  weakly  nonlinear  tidal  systems.  Use  of  processed 
intertidal  hydrodynamics  reduces  computer  disk  space  requirements  by  10-fold  compared 
with  intratidal  hydrodynamics.  Substantial  savings  in  computation  time  for  long-term 
transport  computations  can  also  be  realized. 


Accesion  For 


NTIS  CRA&I  „ 

DTIC  TAB  □ 

Ui. announced  □ 

Justification 


By . 

Dist.ibulion/ 


Availability  Codes 

j  Avail  and  for 
Special 


Unclassified 


security  classification  of  this  page 


C-4nesa-pe<x.}<e  ^stucLrles^  VsJcc4<e.f  rno-sses/ 

to  \a)  ^  M ma~4»' ca  t*  rr>  o  4  s  tc&i 

Sa£}'nily'p  5ea 

O c.ea.n° -frdes>/\j£ix'  ja-i' tons;  T'/’J  a  £  c.u  rren-fs*  cm  mO 


PREFACE 


This  study  was  conducted  as  a  part  of  the  Chesapeake  Bay  water 
quality  model  study,  which  was  contracted  to  US  Army  Engineer  Waterways 
Experiment  Station  (WES)  by  the  US  Environmental  Protection  Agency 
through  the  US  Army  Engineer  District,  Baltimore  (CENAB).  The  Chesa¬ 
peake  Bay  model  study  was  managed  at  WES  by  Mr.  Donald  L.  Robey,  Chief, 
Ecosystem  Research  and  Simulation  Division  (ERSD) ,  Environmental  Labora¬ 
tory  (EL).  The  work  was  monitored  by  Mr.  Larry  Lower  of  CENAB. 

Dr.  John  Harrison  was  Chief  of  the  EL. 

This  study  was  conducted  by  Dr.  Mark  S.  Dortch,  Chief,  Water  Qual¬ 
ity  Modeling  Group  (WQMG) ,  ERSD.  Dr.  Dort-h  prepared  this  report  and 
submitted  it  for  partial  fulfillment  for  the  degree  of  Doctor  of  Philos¬ 
ophy,  Civil  Engineering,  Colorado  State  University  (CSU) ,  Fort  Collins, 
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Thomas  G.  Sanders  for  their  advice  and  recommendations.  Special  thanks 
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port,  and  encouragement.  Dr.  John  M.  Hamrick  of  the  Virginia  Institute 
of  Marine  Science  is  recognized  and  thanked  for  his  valuable  sugges¬ 
tions.  The  author  extends  much  gratitude  to  his  colleague  and  friend, 
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assistance,  recommendations,  and  encouragement  that  he  so  unselfishly 
offered  throughout  this  study. 

Drs.  Billy  H.  Johnson  and  Keu  W.  Kim  of  the  WES  Hydraulics  Labora¬ 
tory  and  Coastal  Engineering  Research  Center,  respectively,  are  thanked 
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and  recognized  for  their  efforts  in  modifying  and  calibrating  the  three- 
dimensional  hydrodynamic  model,  GH3D,  as  part  of  the  Chesapeake  Bay 
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This  report  should  be  cited  as  follows: 

Dortch,  Mark  S.  1990.  "Three-Dimensional,  Lagrangian  Residual 
Transport  Computed  from  an  Intratidal  Hydrodynamic  Model," 

Technical  Report  EL-90- 11,  US  Army  Engineer  Waterways  Experiment 
Station,  Vicksburg,  MS. 
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CHAPTER  1 


INTRODUCTION 


1.1  BACKGROUND 

Considerable  environmental  interest  has  focused  on  pollution  and 
contamination  of  some  of  the  large  estuarine  systems  and  coastal  regions 
of  the  United  States,  as  well  as  the  world.  One  of  the  best  examples  is 
the  concern  for  the  decline  in  the  water  quality  of  Chesapeake  Bay. 

Chesapeake  Bay  is  the  United  States'  largest  estuary,  and  also  one 
of  its  most  productive.  It  supports  important  commercial  and  recrea¬ 
tional  fisheries,  transportation,  industry,  recreation,  and  tourism,  and 
provides  irreplaceable  habitat  for  living  marine  resources  and  wildlife. 
Over  the  past  three  decades,  the  Bay  has  experienced  dramatic  reduc¬ 
tions  in  living  resources,  concurrently  with  decline  in  water  quality 
conditions.  The  US  Environmental  Protection  Agency  (USEPA)  identified 
(USEPA  1983a  and  1983b)  major  contributing  factors  leading  to  the  Bay's 
decline  as  high  concentrations  of  nutrients,  increased  eutrophica¬ 
tion/anoxia,  and  fouling  of  sediments  by  toxic  chemicals. 

Strategies  are  being  sought  to  halt  and  reverse  the  degradation  of 
large,  important  systems  such  as  Chesapeake  Bay.  However,  the  cost  of 
implementing  cleanup  strategies  can  be  formidable.  For  example,  the 
costs  of  cleaning  up  Chesapeake  Bay  pollution  are  estimated  to  be  in  the 
billions  of  dollars.  With  costs  so  high,  every  available  tool  must  be 
used  to  evaluate  a  priori  the  effectiveness  of  proposed  pollution  con¬ 
trol  strategies.  Numerical  water  quality  models,  used  in  conjunction 
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with  monitoring  efforts,  offer  a  relatively  inexpensive  means  of  study¬ 
ing  various  alternatives ,  such  as  nutrient  reduction  goals . 

Monitoring  provides  information  on  past  and  present  state  of  water 
quality.  Monitoring  can  be  used  to  evaluate  past  management  efforts, 
but  it  can  not  be  used  to  estimate  what  future  conditions  might  exist. 
Such  information  can  only  be  obtained  through  technically  sound  model¬ 
ing.  Models  provide  a  flexible,  cost  effective  framework  for  studying 
management  options  and  their  impacts  and  become  the  focal  point  for 
issue  resolution.  Modeling  can  also  help  in  better  understanding  the 
system  and  can  provide  information  for  designing  future  monitoring  pro¬ 
grams  . 

The  physics,  chemistry,  and  biology  of  estuarine/coastal  systems 
are  too  complex  to  base  management  decisions  on  the  results  of  simple 
empirical  or  statistical  models.  Decision  making  can  be  more  soundly 
based  when  information  is  provided  by  mechanistic  simulation  models. 

The  mathematics  of  these  models  are  usually  too  complex  for  analytical 
solutions,  so  mechanistic  simulation  models  are  usually  numerical.  The 
physical,  chemical,  and  biological  processes  to  be  simulated  should  be 
as  well  defined  in  the  model  as  technically  defensible  and  feasible. 

The  development  of  a  numerical  water  quality  model  of  Chesapeake 
Bay  (Dortch  et  al.  1988)  was  initiated  in  1988  to  evaluate  the  future 
effectiveness  of  nutrient  controls  for  improving  water  quality.  This 
model  is  three-dimensional  (3D)  and  time-varying  and  is  coupled  to  a  3D 
hydrodynamic  modal  which  includes  all  the  important  physical  processes 
for  estuaries.  A  bottom  sediment  quality  model  is  coupled  to  the  water 
quality  model  of  the  water  column.  The  sediment  model  simulates  the 
long-term  behavior  of  nutrients  deposited  on  the  bottom  and  their 
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effects  on  the  water  column.  The  Chesapeake  Bay  monitoring  program  was 
tailored  to  provide  information  required  for  the  development  and  cali¬ 
bration  of  the  sediment  quality  model. 

The  research  presented  herein  was  a  critical  component  of  the  over¬ 
all  Chesapeake  Bay  model  development.  This  component  provides  the  in¬ 
terfacing  for  the  water  quality  and  hydrodynamic  models.  Without  the 
proper  interfacing  procedure,  it  would  not  be  feasible  to  apply  the 
water  quality  model  in  a  cost  effective  and  technically  defensible  man¬ 
ner.  Although  advances  in  computer  power  and  speed  have  recently  made 
it  feasible  to  construct  and  apply  time-varying  3D  models,  3D  modeling 
is  still  costly  and  pressing  the  state-of-the-art. 

Three-dimensional,  intratidal  (i.e.  tidally  influenced  or  contains 
tidal  fluctuations)  hydrodynamic  models  typically  have  time  steps  on  the 
order  of  minutes.  For  example,  the  Chesapeake  Bay  hydrodynamic  model 
(HM)  uses  a  time  step  of  five  minutes,  which  is  dictated  by  stability 
requirements.  Ignoring,  for  the  present,  the  stability  requirements  ox 
the  transport  terms  of  the  water  quality  model  (WQM) ,  the  WQM  time  step 
depends  on  the  time  scales  of  the  kinetic  processes,  which  range  on  the 
order  of  hours  to  days.  The  difference  in  the  time  scales  of  the  two 
models  presents  a  problem  since  hydrodynamic  models  are  used  to  drive 
the  transport  terms  of  water  quality  models . 

There  are  basically  two  methods  for  coupling  HM  information  to  the 
WQM;  direct  and  indirect  coupling  (Hall  et  al.  1988).  Direct  coupling 
refers  to  the  use  of  the  same  spatial  grid  and  time  steps  by  both  mod¬ 
els.  Numerical  models  of  lower  dimensionality  carry  sufficiently  low 
computational  burden  to  allow  direct  coupling.  However,  direct  coupling 
may  be  cost  prohibitive  for  models  of  higher  dimensionality  and  greater 
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computational  expense.  In  these  cases,  it  may  be  necessary  to  indirect¬ 
ly  couple  the  two  models.  Indirect  coupling  involves  spatially  and/or 
temporally  averaging  the  HM  output,  storing  this  information,  and  subse¬ 
quently  using  it  as  input  to  drive  the  WQM.  The  temporal  resolution, 
and  possibly  the  spatial  resolution,  necessary  for  hydrodynamic  model 
solutions  may  be  great. c  than  required  for  water  quality  and  can  lead  to 
unacceptable  computational  costs  when  simulating  multiple  water  quality 
constituents  for  long-term  events.  Hydrodynamic  and  water  quality 
models  of  large  tidal  systems  can  be  made  more  tractable  by  indirect 
coupling  of  the  two  models. 

Resolving  long-term  environmental  questions  may  require  simulations 
that  are  impractical  due  to  simulation  costs.  For  example,  the  Chesa¬ 
peake  Bay  model  study  requires  annual  and  multi-year,  even  multi-decade, 
water  quality  simulations  to  properly  evaluate  the  nutrient  reduction 
strategies  (Dortch  et  al.  1988).  The  CPU  requirements  for  an  annual 
water  quality  simulation  of  Chesapeake  Bay  with  intratidal  hydrodynamic 
forcing  are  estimated  to  range  on  the  order  of  hours  on  a  supercomputer. 
The  disk  space  required  to  save  a  year  of  intratidal  hydrodynamic  infor¬ 
mation  is  on  the  order  of  a  billion  bytes  (i.e.  gigabyte)  for  Chesapeake 
Bay.  Developments  in  intertidal  transport  modeling  techniques  can  sig¬ 
nificantly  reduce  these  requirements . 

Averaging  tidally  varying  HM  output  over  periods  on  the  order  of 
the  tidal  period,  or  longer,  produces  intertidal  (residual)  currents 
that  have  considerably  less  magnitude  than  intratidal  currents  (i.e. 
currents  averaged  over  time  periods  less  than  a  tidal  period) .  For 
example,  the  tidally  averaged  current  of  a  truly  periodic  flow  is  zero. 
The  use  of  residual  currents  can  significantly  reduce  the  stability 
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requirements  for  explicitly  computed  advective  flow  terms  of  the  water 
quality  transport  equation,  thus,  allowing  larger  WQM  time  steps.  Large 
amounts  of  disk  space  can  be  easily  consumed  when  storing  time -varying, 
3D  velocities.  Additionally,  reading  in  large  amounts  of  HM  output  dur¬ 
ing  WQM  execution  can  significantly  slow  down  computation  speed.  The 
use  of  residual  (intertidal)  velocities,  as  opposed  to  intratidal  veloc¬ 
ities,  reduces  these  requirements  by  about  an  order  of  magnitude  since 
intertidal  information  updates  are  on  the  order  of  12.4  hours  or  more, 
whereas  intratidal  information  updates  are  on  the  order  of  one  to  two 
hours.  Therefore,  the  use  of  intertidal  transport  for  the  WQM  can  sig¬ 
nificantly  reduce  computational  expense,  making  multiple  constituent, 

3D,  long-term  water  quality  simulation  costs  more  reasonable.  The  pur¬ 
pose  of  the  research  presented  herein  is  to  develop  residual  (inter¬ 
tidal)  transport  modeling  to  reduce  computational  and  disk  storage  re¬ 
quirements,  while  retaining  tidal  influences. 

1 . 2  PROBLEM 

There  are  problems  with  computing  the  residual  currents,  from  the 
basic  intratidal  HM  information,  in  a  manner  that  preserves  the  correct 
transport  characteristics.  Simply  averaging,  over  one  or  more  tidal 
periods,  the  HM  velocities  at  each  point  produces  Eulerian  residual  cur¬ 
rents.  The  mass  passing  a  fixed  point  may  not  depend  solely  on  the  mean 
velocity  at  that  point,  but  it  may  depend  on  other  properties  of  the 
flow  field,  such  as  the  interactions  with  neighboring  velocities.  The 
correct  residual  velocities  for  mass  transport  are  of  a  Lagrangian 
nature  (Longuet-Higgins  1969,  Cheng  and  Casulli  1982,  Feng  et  al.  1986a 
and  1986b,  and  Orbi  and  Salomon  1988),  i.e.  the  net  displacement  of 
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marked  water  parcels  divided  by  the  elapsed  time,  or  the  time  average  of 
the  instantaneous  velocity  of  a  particle.  Lagrangian  residual  currents, 
which  can  be  quite  different  from  Euler ian  residual  currents,  are  a 
result  of  interactions  of  system  forcing  (inflows,  tides,  wind,  and 
earth's  rotation)  with  system  characteristics  (geometry,  bathymetry,  and 
density  gradients) . 

Although  there  is  a  recognized  need  for  the  use  of  Lagrangian  re¬ 
sidual  currents  in  intertidal  transport  modeling,  there  are  no  known 
examples  of  computing  3D  Lagrangian  residual  circulation  from  an  in- 
tratidal  HM  for  use  in  a  time -varying,  intertidal  transport  model. 
Therefore,  practical  information  is  not  available  for  implementing  3D 
Lagrangian  residual  transport. 

1.3  RESEARCH  OBJECTIVES 

The  basic  hypothes-'s  of  this  research  is  that  tide-induced  residual 
currents  can  not  be  ignored  in  residual  transport  modeling  of  Chesapeake 
Bay.  The  goal  of  this  research  is  to  develop  a  method  for  computing  3D 
Lagrangian  residual  currents  from  an  intratidal  hydrodynamic  model  for 
use  in  an  intertidal  transport  model.  The  following  objectives  have 
been  established  to  accomplish  this  goal. 

1.  Indirectly  couple  (interface)  the  HM  and  WQM  such  that  the 
transport  characteristics  of  the  HM  are  preserved  in  the  WQM  for 
intratidal  averaging  periods  (i.e.  one-  or  two-hour  averaging  in- 
tsrvsls  or  loss) » 

2.  Develop  and  implement  the  Lagrangian  residual  computations  in  a 
way  that  will  ensure  mass  conservation  when  applied  to  the 
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intertidal  transport  equation.  The  formulation  must  be  compatible 
with  the  hydrodynamic  and  transport  modeling  frameworks. 

3.  Verify  the  computational  procedure  by  comparing  numerical  re¬ 
sults  with  the  two-dimensional,  analytical  results  of  Iannello 
(1977). 

4.  Apply  the  procedure  to  Chesapeake  Bay  and  evaluate  intertidal 
transport  through  comparisons  to  observed  salinity  data  and  salini¬ 
ty  computed  by  the  HM. 

5.  Investigate  the  characteristics  of  3D  Lagrangian  residual  cur¬ 
rents  . 

The  basic  theory  and  formulations  are  presented  in  Chapter  2;  the 
computational  procedures  are  developed  in  Chapter  3;  and  the  methodology 
is  applied  and  evaluated  in  Chapter  4. 


CHAPTER  2 


LITERATURE  REVIEW 

This  chapter  presents  the  theory  and  basic  formulations  for  comput¬ 
ing  residual  currents.  The  term  residual  currents  refers  to  filtering 
or  averaging  out  the  intratidal  fluctuations.  In  the  first  section,  the 
tidally  averaged  (intertidal)  transport  equation  is  obtained,  and  the 
concept  of  tidal  dispersion  used  with  Eulerian  residual  circulation  is 
discussed.  The  concept  of  Lagrangian  residual  circulation  is  then  de¬ 
scribed,  and  the  Lagrangian  intertidal  transport  equation  is  presented. 
Approximations  for  Lagrangian  residuals  are  discussed.  The  first-order 
approximation  is  the  sum  of  the  Eulerian  residual  and  the  Stokes'  drift. 
Stokes'  drift  is  defined,  and  both  the  original  and  mass  conserving 
formulations  are  presented  in  Section  2.3.  Previous  studies  of  Lagran¬ 
gian  residual  circulation  are  discussed  in  Section  2.4.  The  chapter  is 
summarized  in  the  last  section. 

2.1  INTERTIDAL  TRANSPORT  AND  EULERIAN  RESIDUALS 

The  basis  for  mechanistic  water  quality  models  (and  other  types  of 
transport  models)  is  the  conservation  of  mass,  or  mass  transport  equa- 
tior .  The  3D  mass  transport  equation  in  car-tesian  coordinates  and  ten¬ 
sor  notation  (Yih  1977)  is 

*  5  SOURCES/SINKS  (2.1) 
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where 

G  -  constituent  concentration 
DiJ  -  turbulent  eddy  diffusivity  Coefficient 

t  -  time 

UA  -  flow  velocity  in  direction  i 

XA  -  coordinate  in  direction  i 

The  first  term  in  Equation  2.1  is  the  time  rate  of  change  of  mass  con¬ 
centration  within  a  fluid  element.  The  second  term  is  the  advection  of 
mass  per  unit  volume  resulting  from  flow  into  and  out  of  the  fluid  ele¬ 
ment.  The  third  term  (first  term  on  the  right  side  of  Equation  2.1)  is 
the  diffusion  of  mass  per  unit  volume  across  the  boundaries  of  the  fluid 
element.  This  term  usually  includes  molecular  and  turbulent  diffusion, 
and  it  may  also  include  shear  dispersion  for  models  that  average 
spatially  in  one  or  more  dimensions.  The  last  term  represents  the  rate 
at  which  mass  per  unit  volume  is  added  to  (sources)  and/or  taken  from 
(sinks)  the  fluid  element  by  various  internal  transfers  and  transforma¬ 
tions.  Equation  2.1  is  solved  for  each  water  quality  state  variable 
(dependent  variable)  over  a  specified  spatial  and  temporal  domain 
(independent  variables).  The  velocities  and  diffusivities  are  a  result 
of  the  hydrodynamics  of  the  system  and  play  an  important  role  in  deter¬ 
mining  the  transport  of  salinity,  sediment,  nutrients,  and  other  dis¬ 
solved  or  suspended  matter.  The  hydrodynamic  information  used  to  drive 
the  transport  model  is  usually  obtained  from  a  hydrodynamic  model. 

Water  quality  model  studies  of  tidal  systems  commonly  use  inter¬ 
tidal  (residual)  mass  transport.  This  means  that  the  mass  concentra¬ 
tions  are  either  steady-state  or  slowly  time-varying  for  periods  on  the 
order  of  a  tidal  cycle.  The  hydrodynamics  used  to  drive  the  WQM  are 
averaged  over  a  tidal  cycle.  The  primary  advantage  of  this  approach  is 
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that  the  intertidal  transport  filters  out  tidal  fluctuations  and  focuses 
on  long-term  variations  resulting  from  mean  flow  considerations.  Nihoul 
and  Ronday  (1975)  and  others  have  pointed  out  that  residual  currents, 
rather  than  tidal  currents,  determine  the  overall  ecological  balance  or 
the  long-term  movement  of  water  quality  constituents  and  pollutants  in 
tidal  systems. 

The  tidally  averaged  transport  equation  is  obtained  by  first  decom¬ 
posing  (Officer  1976)  instantaneous  variables  into  tidally  averaged  and 
tidally  varying  components,  i.e. 

*-?  +  *'  (2.2) 

where 

1  t  +T 

-  4  /  0  <f>  dt  -  tidally  averaged  variable 

i  t0 

V  -  o 

and  T  is  the  averaging  period  (e.g.  tidal  period  or  longer).  Implement¬ 
ing  Equation  2.2  for  U  and  C  in  Equation  2.1,  ignoring  sources  and 
sinks,  and  averaging  over  a  tidal  cycle  (recognizing  that  the  averages 
of  all  cross  product  terms  involving  mean  and  tidally  fluctuating  com¬ 
ponents  are  zero  and  neglecting,  for  advection  dominated  systems,  tidal¬ 
ly  fluctuating  correlations  of  diffusivity  and  concentration)  results  in 
the  tidally  averaged  transport  equation 

ac  .  3(0,0)  .  3(uTT)  5(DiJ  "ax;] 
at  T  ax,  ‘r  ax,  “  ax, 

The  mean  velocities  in  the  second  term  of  Equation  2.3,  referred  to 
as  Eulerian  residual  velocities  (Officer  1976),  Ug,  are  obtained  by 
averaging  the  velocity  at  fixed  points  over  one  or  more  tidal  cycles. 
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The  third  term  in  Equation  2.3,  which  is  the  non-zero  correlation  be¬ 
tween  the  tidal  velocity  and  concentration  fluctuations  as  a  result  of 
time  averaging  of  the  advective  terms,  has  been  referred  to  as  tidal 
dispersion  (Officer  1976,  Fischer  et  al.  1979).  Assuming  the  Fickian 
form  of  diffusion  for  mean  concentration,  the  third  term  can  be  rewrit¬ 
ten  as 


-U'C' 


(2.4) 


where  Dr  is  the  tidal  dispersion  coefficient.  It  has  been  common  prac¬ 
tice  in  estuarine  transport  modeling  to  add  tidal  dispersion  to  the 
diffusion/dispersion  term  on  the  right  side  of  Equation  2.3,  resulting 
in  a  total  dispersion  term  that  includes  turbulent  diffusion,  shear  dis¬ 
persion  (for  spatially  averaged  models) ,  and  tidal  dispersion  (Fischer 
et  al.  1979).  Therefore,  Equation  2.3,  expressed  in  terms  of  mean  vari¬ 
ables  only,  becomes  the  Eulerian  residual  transport  equation, 


ac  ,  8<V>  3(Dl“  8*J 

at  aXi  "  dxL 


(2.5) 


There  are  numerous  examples  of  the  use  of  Eulerian  residual  cir¬ 
culation  for  driving  water  quality  transport  models;  a  recent  example  is 
the  steady-state  Chesapeake  Bay  water  quality  model  (HydroQual  1987  and 
Fitzpatrick  et  al.  1988).  The  problem  with  this  approach  is  that  an 
advective  process  is  lumped  into  Fickian  diffusion  terms,  resulting  in 
an  unrealistic  representation.  There  has  been  large  variability  in  the 
range  of  observed  tidal  dispersion  coefficients  since  it  incorporates 
tidal  fluctuations  that  can  vary  widely  in  space  and  time. 
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Additionally,  the  turbulent- diffusion  and  shear  dispersion  become  rela¬ 
tively  unimportant  when  introducing  tidal  dispersion  (Dyer  1973,  Fischer 
1976). 

2.2  LAGRANGIAN  RESIDUALS 

Researchers  (Longuet-Higgins  1969,  Zimmerman  1979,  Awaji  1982, 

Cheng  and  Casulli  1982,  Feng  et  al.  1986a  and  1986b,  and  Orbi  and 
Salomon  1988)  have  recognized  the  need  to  use  Lagrangian  residual  cur¬ 
rents  rather  than  Eulerian  residual  currents  to  properly  describe  the 
origin  of  water  masses.  Lagrangian  residual  currents  are  related  to 
Lagrangian  mean  velocities ,  which  are  the  average  velocities  of  marked 
water  parcels  tracked  over  one  or  more  tidal  cycles  (Feng  1987).  The 
Lagrangian  mean  velocity  is  also  described  as  the  net  displacement  of  a 
marked  particle  over  one  or  more  tidal  cycles  divided  by  the  displace¬ 
ment  time.  Feng  (1987)  points  out  that  the  Lagrangian  residual  velocity 
can  be  defined  as  the  Lagrangian  mean  velocity  and  can  be  used  as  an 
Eulerian  field  variable  in  the  mass  transport  equation  if  the  Lagrangian 
mean  velocities  satisfy  continuity.  Such  treatment  eliminates  the  need 
to  include  the  tidal  dispersion  effect  in  the  diffusion/dispersion  terms 
of  the  transport  equation,  at  least  for  weakly  nonlinear  systems  (Feng 
et  al.  1986b). 

The  intertidal  (residual)  mass  transport  equation  has  been  derived 
for  two-dimensional,  depth  integrated  flow  and  three-dimensional  flow  by 
Feng  et  al.  (1986a  and  1986b)  and  Feng  (1987),  respectively.  A  small 
parameter  perturbation  technique  (Van  Dyke  1964)  and  tidal  averaging 
were  used  to  develop  the  solution  from  the  intratidal  nondimens ional 
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transport  equation.  The  small  parameter,  «,  is  a  measure  of  the  non¬ 
linearity  of  the  system, 


where 


(2.6) 


and  where  the  characteristic  values  are  denoted  as:  tidal  amplitude,  J’c; 
water  depth,  hc;  tidal  excursion,  lc;  and  basin  horizontal  length  scale, 
Lc.  The  solution  is  valid  for  weakly  nonlinear  systems,  i.e.  small  k  or 
k  <  1.0.  With  the  expansion  solution  carried  to  order  k?  ,  the  steady  3D 
residual  transport  equation  derived  by  Feng  (1987)  is  stated  as 


U, 


LM 


(2.7) 


where 

ULM  -  the  3D,  first-order  Lagrangian  residual  velocity  vector 
V  -  gradient  operator 

C  -  the  intertidal,  long-term  average  concentration 
Dz  -  tidally  averaged  vertical  eddy  diffusivity 
z  -  vertical  coordinate 


and  bold  characters  represent  vector  quantities.  Equation  2.7  has  the 
form  of  the  Eulerian  residual  transport  equation  (Equation  2.5),  the 
primary  difference  being  that  Lagrangian  residual  velocities  are  used  as 
Eulerian  field  variables  rather  than  Eulerian  residual  velocities.  It 
should  also  be  noted  that  Equation  2.7  does  not  contain  the  tidal  dis¬ 
persion  terms  presented  in  Equation  2.5.  The  Lagrangian  residual  veloc¬ 
ities  have  included  the  effect  of  the  tidal  fluctuations,  at  least  to  a 
first-order  approximation.  The  first-order  Lagrangian  residual  velocity 
is  the  sum  of  the  Eulerian  residual  velocity  and  the  Stokes'  drift, 
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which  will  be  discussed  in  Section  2.3.  Horizontal  diffusion,  which  is 
not  included  in  Equation  2.7,  only  provides  higher  order  accuracy  to  the 
advection  dominated  transport  (Feng  1986b  and  1987). 

A  result  similar  to  Equation  2.7  was  obtained  earlier  by  Andrews 
and  McIntyre  (1978).  Through  Eulerian-Lagrangian  transformation  theory, 
Andrews  and  McIntyre  developed  an  exact  theory  for  generalized  Lagran- 
gian-mean  flow  subject  to  finite -amplitude  disturbances.  The  resulting 
"generalized  Lagrangian-mean"  (GLM)  operator  describes  wave  and  mean 
flow  interactions  with  equations  in  Euler ian  form.  The  GLM  operator  for 
the  material  derivative  is  stated  as 

dL  ”  it  +  '  7  (2.8) 

where  UL  is  the  Lagrangian  residual  velocity  vector.  Equation  2.8  was 
obtained  by  requiring  that  the  mean  of  the  disturbance-associated  par¬ 
ticle  displacement  field  is  zero,  i.e. 

-  0.  (2.9) 

If  £  is  the  displacement  associated  with  the  fluctuating  tidal  veloci¬ 
ties,  U' ,  then  U'  -  0,  which  is  the  case  for  instantaneous  velocities 
decomposed  into  tidally  averaged  and  fluctuating  components.  Andrews 
and  McIntyre  (1978)  indicate  that  the  difference  between  the  Lagrangian 
mean  description  (Equation  2.8)  and  the  Eulerian  mean  description  (i.e. 
left  side  of  Equation  2.5)  is  accounted  for  by  the  Stokes'  correction. 

Hamrick  (1987)  developed  a  3D  tidally  averaged  mass  transport  equa¬ 
tion  for  a  vertically  stretched  and  horizontally  curvilinear  boundary- 
fitted  coordinate  system  using  small  parameter  perturbation  techniques. 
These  results  are  useful  in  this  research  since  the  hydrodynamic  model 
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that  is  used  is  based  on  boundary- fitted  coordinates.  Hamrick's  solu¬ 
tions,  which  were  carried  to  the  same  order  as  Feng's  (1987)  results, 
also  confirm  that  the  intertidal  transport  equation  is  driven  by  Lagran- 
gian  residual  velocities,  which  are  equal  to  the  sum  of  the  Eulerian 
residuals  and  the  Stokes'  drift  at  the  first-order  of  approximation. 
Additionally,  Hamrick's  solutions  resulted  in  3D  intertidal  diffusivi- 
ties  that  were  simply  intertidal  means  (i.e.  time  averages  for  one  or 
more  tidal  cycles  or  low  pass  filters)  of  the  instantaneous  turbulent 
diffusivities .  The  implications  of  boundary-fitted  coordinates  will  be 
discussed  in  Chapter  3.  For  now,  the  general  time -varying,  3D  Lagran- 
gian  residual  transport  equation  in  tensor  notation  for  cartesian  coor¬ 
dinates  is  stated  as 


Earlier,  Zimmerman  (1979)  provided  a  substantial  improvement  in 
understanding  Lagrangian  residual  currents.  He  showed  through  an  Euler- 
Lagrangian  transformation  that  the  Stokes'  drift  was  only  a  first-order 
approximation  of  the  Lagrangian  residual  current.  The  perturbation  ana¬ 
lysis  by  Feng  et  al.  (1986a)  confirmed  that  the  first-order  Lagrangian 
residual  current  is  the  sum  of  the  Eulerian  residual  current  and  the 
Stokes'  drift.  The  second-order  solution  of  Feng  et  al.  (1986a)  shows 
that  the  Lagrangian  residual  includes  a  Lagrangian  drift  term  that  is  a 
periodic  function  of  the  tidal  phase, 

ut  "  ue  +  us  +  *  ULD  (2 . 13.) 


where 
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UE  -  Eulerian  residual  velocity 
Us  -  Stokes'  drift  velocity 
ULD  -  Lagrangian  drift  velocity 

k  -  measure  of  system  nonlinearity,  defined  previously 
The  Lagrangian  drift  velocity  is  induced  by  nonlinear  interactions  be¬ 
tween  tidal  currents  and  the  first-order  residual  currents.  Conse¬ 
quently,  the  second-order  Lagrangian  residual  circulation  depends  on  the 
particle  release  time  and  is  tidal  phase  dependent. 

Figure  2.1  demonstrates  the  concepts  of  residual  currents.  Con¬ 
sider  a  tidal  basin  where  particles  can  be  released  and  tracked  over  a 
tidal  cycle.  A  particle  is  released  at  time  tQ  and  is  tracked  through  a 
complete  tidal  cycle.  The  terminus  of  this  particle  is  displaced  from 
the  release  point  indicating  a  net  drift.  If  other  particles,  which  are 
released  at  other  times  within  the  tidal  cycle  (e.g.  one  hour  inter¬ 
vals),  arrive  at  the  same  end  point,  the  distance  between  the  release 
point  and  the  termini  represents  the  net  displacement  associated  with 
the  Eulerian  residual  and  Stokes'  drift.  The  Lagrangian  drift  is  negli¬ 
gible  in  this  case.  If  particles  released  at  different  times  within  the 
tidal  cycle  (e.g.  tQ  +  1.0  hour)  have  different  termini,  then  Lagran¬ 
gian  drift  is  evident.  The  distance  between  the  termini  of  the  parti¬ 
cles  and  the  centroid  of  their  termini  represent  Lagrangian  drift. 

It  stands  to  reason  that  if  a  system  is  very  weakly  nonlinear  (i.e. 
very  small  k) ,  then  second-order  tidal  phase  effects  are  negligible,  and 
the  first-order  approximation  for  Lagrangian  residuals  is  sufficient. 

The  first-order  approximation  may  not  be  sufficient  for  systems  with 
considerable  nonlinearity.  The  English  Channel  is  a  good  example  of  a 
tidal  system  with  considerable  nonlinearity  (Orbi  and  Salomon  1988) , 
where  k  approaches  0.25.  The  calculation  of  tracer  trajectories  by 
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Cheng  (1983)  with  a  2D  depth -averaged  model  of  South  San  Francisco  Bay 
(k  between  0.2  and  0.5)  showed  that  the  computed  Lagrangian  current  de¬ 
pended  on  the  particle  release  time.  In  these  cases,  second-order  tidal 
phasing  effects  may  need  to  be  taken  into  account,  along  with  a  rela¬ 
tively  fine  spatial  grid  resolution  (Orbi  and  Salomon  1988  and  Cheng 
1983). 

A  mathematical  relationship  for  computing  Ujjj  was  determined  by 
Feng  et  al.  (1986a)  for  an  analytical  test  case  of  an  M2  tide.  An  M2 
tide  is  a  12.42  hour  period  harmonic  tidal  wave  constituent  associated 
with  the  principal  lunar  component.  For  real  tidal  embayments,  there 
are  no  analytical  solutions.  Particle  tracking  techniques  for  succes¬ 
sively  released  particles  (e.g.  released  at  one  hour  intervals)  could  be 
used  to  determine  tidal  phase  dependent  Lagrangian  displacements,  thus 
obtaining  tidally  varying  Lagrangian  velocities  (Cheng  1983  and  Orbi  and 
Salomon  1988).  However,  such  intratidal  Lagrangian  information  is  not 
consistent  with  the  interest  here,  i.e.  to  use  intertidal  residual  cur¬ 
rents  for  long-term  transport.  It  does  seem  feasible  to  compute  the 
center  of  mass  of  particles  released  throughout  a  tidal  cycle  to  obtain 
the  tidally  averaged  Lagrangian  displacements  (thus  the  Lagrangian  re¬ 
siduals).  The  spread  of  the  tidally  averaged  trajectories  of  the  parti¬ 
cles  could  provide  an  estimate  of  the  tidal  phase  induced  dispersion 
that  arises  from  tidally  averaging  a  nonlinear  system  (Cheng  1983) . 
Computation  of  Lagrangian  residual  currents  that  include  the  second-or¬ 
der  tidal  phase  effects  would  be  a  logical  extension  in  future  research, 
but  it  is  beyond  the  scope  of  this  research. 

In  this  work,  the  Lagrangian  residual  circulation  will  be  computed 
using  the  first-order  approximation  (Equation  2.11  without  the 
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terra).  The  theory  (Feng  et  al.  1986a  and  1986b  and  Hamrick  1987)  indi¬ 
cates  that  this  approximation  should  be  sufficiently  accurate  for  the 
intended  use  (i.e.  weakly  nonlinear  tidal  embayments,  such  as  Chesapeake 
Bay) .  It  is  reasonable  to  assume  that  Chesapeake  Bay  is  weakly  non¬ 
linear  since  the  maximum  tidal  amplitude  of  approximately  0.4  m  and  mean 
depth  of  about  8.0  m  (Fisher  1986)  yield  a  k  of  0.05  or  less. 

2.3  STOKES'  DRIFT 

Stokes'  drift  is  the  correction  velocity  at  a  fixed  point  that  is 
added  to  the  Eulerian  residual  to  produce  the  first-order  estimate  of 
the  Lagrangian  residual.  A  formulation  for  Stokes'  drift  is  necessary 
to  determine  the  Lagrangian  residual  velocities  which  are  used  within 
the  Lagrangian  residual  transport  equation  (Equation  2.10). 

An  original  formulation  for  Stokes'  drift  is  presented  first  to 
examine  the  physical  meaning  of  Stokes'  drift.  However,  the  original 
formulation  does  not  guarantee  preservation  of  continuity  (i.e.  flow  and 
volume  conservation).  A  mass  conserving  formulation  follows  the  ori¬ 
ginal  formulation. 

2.3.1  Original  Formulation 

The  first  use  of  the  concept  of  Lagrangian  residual  transport  was 
by  Longuet-Higgins  (1969)  who  derived  the  first-order  of  approximation 
for  the  Lagrangian  residual,  which  is  equal  to  the  sum  of  the  Eulerian 
residual  velocity  and  the  Stokes'  drift.  Longuet-Higgins  started  with  a 
first-order  Taylor  series  expansion  for  the  Eulerian  velocity  field  to 
describe  the  velocity  of  a  particle  in  an  oscillating  flow, 


U(X,t)  -  U(X0,t)  +  AX  •  V  U(X0,t) 


(2.12) 
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where , 

U(X,t)  -  particle  velocity  at  new  location,  X 
X  -  X0  +  AX  «  particle  location  at  time  t 
U(X0,t)  -  particle  velocity  at  old  location,  X0 

If  AX  is  small  compared  to  the  local  length-scale  of  the  velocity  field, 
then  AX  may  be  approximated  as 

AX  -  J*  U(XQ,t)  dt  (2.13) 

to 

Substituting  Equation  2.13  into  2.12  and  taking  mean  values  over  one  or 
more  tidal  cycles  results  in 

U(X,t)  -  U(X0,t)  +  f  U(X0,t)  dt  •  V  U(X0,t)  (2.14) 

where  the  overbars  represent  time  averaging.  Longet-Higgins  referred  to 
the  left-hand  side  as  the  mass  transport  velocity,  which  is  the  sum  of 
the  Eulerian  residual  velocity  and  the  Stokes'  drift;  thus,  the  Stokes' 
drift  is 

Us  -  J  U  dt  •  V  U  (2. 15) 

The  Stokes'  drift  can  be  thought  of  as  the  residual  current  resulting 
from  the  time  average  of  the  spatial  variability  of  the  Eulerian  veloc¬ 
ity  field.  Stokes'  drift  is  induced  from  the  nonlinear  interaction  of 
the  tidal  currents  (Feng  et  al.  1986a). 

The  velocities  in  Equation  2.15  are  the  total  instantaneous  veloci¬ 
ties  for  an  oscillating  flow.  Flow  fields  containing  mean  flow  compo¬ 
nents  can  be  decomposed  into  the  mean  and  tidal  fluctuating  components 
(see  Equation  2.2).  According  to  the  theory  of  Andrews  and  McIntyre 
(1978),  the  velocities  in  Equation  2.15  should  be  the  tidally 
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fluctuating  components,  U' .  The  perturbation  analyses  by  Feng  et  al. 
(1986a  &  b)  and  Hamrick  (1987)  also  indicate  that  periodic  components 
are  appropriate  for  computing  the  Stokes'  drift.  Following  the  theory 
and  results  of  the  latter  research,  the  generalized  Stokes'  drift  for¬ 
mulation  is  stated  as 

Us  -  J  U'  dt  •  V  U'_  (2.16) 

2.3.2  Mass  Conserving  Formulation 

The  above  formulation  for  Stokes'  drift  (Equation  2.16)  may  not 
guarantee  mass  conservation  when  implemented  in  a  numerical  calculation. 
An  alternate  form  of  Equation  2.16  can  be  obtained  which  will  guarantee 
mass  conservation  (Hamrick  1987). 

Following  from  Longuet-Higgins '  (1969)  analysis  of  a  periodic  cur¬ 
rent,  Equation  2.16  can  be  written  in  the  form 

Us  -  7  x  B  -  curl  B  (2.17) 

where  the  components  of  B  are  defined  as 

Bx  -  v'  J  w'  dt 
By  -  w'  J  u'  dt 

Bz  -  “J  v'  dt  (2.18) 

Equations  2.17  and  2.18  can  be  derived  from  Equation  2.16  through  use  of 
the  continuity  equation  and  the  fact  that  if  A  and  B  are  any  two  peri¬ 
odic  quantities  with  zero  mean  (e.g.  U'  -0),  then 

* 


A  /  B  dt  +  B  /  A  dt  -  0. 


(2.19) 
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Equation  2.17  ensures  mass  conservation  (e.g.  V  •  Ug  -  0)  since  the  di¬ 
vergence  of  a  curl  is  zero  (Sears  1970) .  If  the  Eulerian  residuals  are 
obtained  from  a  hydrodynamic  model  that  is  based  on  conservation  of  mass 
and  momentum,  then  V  •  Ug  -  0.  Therefore,  with  conservative  Eulerian 
residuals  and  Stokes'  drifts,  the  first-order  Lagrangian  residuals  will 
also  conserve  mass. 


2.4  PREVIOUS  STUDIES  OF  LAGRANGIAN  RESIDUALS 


This  section  reviews  previous  studies  of  Lagrangian  residual  cir¬ 
culation.  Although  a  number  of  researchers  have  investigated  Lagrangian 
residual  currents,  the  work  to  date  on  application  of  Lagrangian  resid¬ 
ual  circulation  in  transport  modeling  has  been  very  limited. 

Tee  (1976)  was  probably  the  first  to  use  a  numerical  model  to  com¬ 
pute  residual  currents.  He  processed  Eulerian  residuals  from  a  2D 
depth- averaged  nonlinear  hydrodynamic  model.  He  also  attempted  to  com¬ 
pute  the  Stokes'  drift;  however,  his  Stokes'  drift  was  the  difference  in 
the  Eulerian  mean  volumetric  transport  divided  by  the  mean  water  depth 
and  the  Eulerian  residual  velocities.  The  Stokes'  drift  given  by  Tee 
(1976)  is  only  valid  for  ID  flows  (Feng  et  al.  1986a).  Tee  did  not  use 
the  computed  residual  currents  for  transport  simulations. 

Cheng  and  Casulli  (1982)  provided  considerably  improved  insight 
into  the  nature  of  Eulerian  and  Lagrangian  residual  currents.  They  ap¬ 
plied  a  2D  depth- integrated  model  to  South  San  Francisco  Day.  The  model 


of  the  bay  was  driven  to  a  dynamic  steady-state  with  art  M2  tide.  Euler¬ 


ian  residual  currents,  which  were  computed  from  tidally  averaged  hydro- 
dynamic  output,  were  found  to  be  quite  different  from  Lagrangian  resid¬ 
uals,  which  were  obtained  through  particle  tracking  (e.g.  particle 
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displacement  over  a  tidal  cycle  divided  by  the  tidal  period) .  Cheng  and 
Casulli  also  found  that  in  some  areas  the  Lagrangian  residuals  depended 
on  the  particle  release  time  or  the  phase  of  the  tide.  Tidal  phase  de¬ 
pendency  is  not  surprising  in  South  San  Francisco  Bay  since  major  por¬ 
tions  are  relatively  shallow  (i.e.  less  than  2.0  m) ,  and  the  tidal  am¬ 
plitude  is  rather  large  (on  the  order  of  1.0  m) ;  thus,  areas  of  South 
San  Francisco  have  strong  topographic  influence  with  considerable  non¬ 
linearity.  Cheng  and  Casulli  (1982)  did  not  use  their  computed  residual 
currents  to  drive  a  transport  model. 

Awaji  (1982)  also  used  particle  tracking  to  study  Lagrangian  move¬ 
ment  through  a  coastal  strait  with  and  without  the  effect  of  turbulence 
generated  by  a  Markov-chain  random  walk  procedure.  He  used  a  2D  depth- 
integrated  hydrodynamic  model  driven  with  an  M2  tide  to  develop  dynamic 
steady-state  tidal  currents  within  an  outer  and  inner  bay  co:  'ected  by  a 
narrow  strait.  Particles  released  throughout  the  grid  were  tracked  over 
three  tidal  cycles  with  and  without  random  turbulent  velocities.  Al¬ 
though  the  imposed  turbulence  created  final  trajectories  that  were  dif¬ 
ferent  from  those  without  turbulence,  the  study  did  not  provide  clear 
insight  into  the  relationship  of  tidal  exchange,  mixing,  and  residual 
currents.  Awaji  (1982)  did  not  compute  residual  currents  in  his  study. 

Cheng  (1983)  recognized  that  the  time  scale  for  ecological  pro¬ 
cesses  is  much  longer  than  the  tidal  period,  and  it  is  impractical  to 
formulate  ecological  models  on  the  same  time  scale  as  tidal  circulation 
models.  He  also  recognized  the  Lagrangian  nature  of  transport  phenom¬ 
ena.  Using  the  dynamic  steady-state  Eulerian  flow  field  generated  with 
a  2D  depth- averaged  hydrodynamic  model  of  South  San  Francisco  Bay,  Cheng 
(1983)  extended  the  results  of  Cheng  and  Casulli  (1982)  by  computing 
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bay-wide  Lagrangian  residual  circulation  by  tracking,  throughout  a  tidal 
cycle,  particles  released  at  one  hour  intervals  over  a  semi-diurnal  tide 
of  12  hour  period.  At  most  of  the  release  points,  the  Lagrangian  resid¬ 
ual  currents  depended  on  the  release  time,  or  the  phase  of  the  tide, 
which  is  not  surprising  considering  the  rather  large  degree  of  nonline¬ 
arity  of  this  system  as  mentioned  above.  The  spread  of  the  Lagrangian 
residual  vectors  indicate  a  mechanism  of  tidal  current  phase  induced, 
mixing.  Cheng  (1983)  did  not  attempt  to  use  these  results  or  methods  to 
drive  a  transport  model. 

Cheng  et.  al.  (1984)  used  an  Eulerian-Lagrangian  Method  (ELM)  to. 
transport  salt  in  South  San  Francisco  Bay.  The  ELM  uses  Lagrangian  par¬ 
ticle  tracking  and  interpolation  to  the  fixed  Eulerian  grid  points  for 
the  advective  transport.  Usually  all  other  processes  (diffusion  and 
transformation),  are  computed  on  the  Eulerian  reference  frame.  Thus,  ELM 
has  two  basic  components,  particle  tracking  and  interpolation  onto  the 
Eulerian  grid.  Cheng  el  al.  (1984)  used  the  2D  dynamic  steady-state 
flow  field  discussed  in  the  previous  paragraph  to  drive  salinity  trans¬ 
port  using  the  ELM.  The  salinity  transport  model  was  intratidal,  thus-, 
residual  currents  were  not  computed  nor  used  for  transport.  The  paper 
focused  on  higher  order  interpolation  techniques  to  reduce  artificial 
numerical  diffusion;  mass  conservation  properties  of  the  ELM  were  not 
discussed. 

The  advantages  of  the  ELM  are:  particle,  tracking  can  provide  a 
direct  means  of  computing  Lagrangian  residual  circulation;  ELM  can  re¬ 
duce  numerical  stability  requirements  since  the  Eulerian  advection 
scheme  is  removed  (Baptista  et  al.  1984);  and  if  high  order  grid  inter¬ 
polations  are  used,  ELM  can  reduce  numerical  dampening  associated  with 
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low  order  (i.e.  upwind  differencing)  Eulerian  advection  schemes  (Cheng 
et  al  1984).  The  disadvantages  of  ELM  are:  computational  techniques 
based  on  the  Lagrangian  viewpoint  are  not  as  well  developed  and  as  ad¬ 
vanced  as  Eulerian  methods  (Cheng  1983),  and  3-D  particle  tracking  for 
long-term  simulations  would  pose  a  computational  challenge;  ELM  poses 
problems  near  boundaries  when  the  particle  trajectory  extends  outside  of 
the  flow  domain;  and  ELM  can  not  guarantee  mass  conservation  (Benque  et 
al.  1982),  whereas  pure  Eulerian  methods  can.  The  lack  of  mass  conser¬ 
vation  may  not  be  a  problem  for  short-term  simulations,  but  it  could  be 
a  serious  limitation  for  long-term  water  quality  simulations.  For  this 
reason,  this  research  focuses  on  Eulerian  methods  for  computing  and  ac¬ 
complishing  residual  transport.  Also,  it  should  be  noted  that  the  sta¬ 
bility  requirements  imposed  by  explicit  Eulerian  advection  schemes  are 
greatly  reduced  for  intertidal  residual  currents  because  of  the  small 
magnitude  of  these  currents. 

The  work  of  Feng  et  al.  (1986a  and  1986b)  contributed  significantly 
to  understanding  tide-induced  Lagrangian  residual  currents  and  their  use 
in  intertidal  transport.  Using  perturbation  techniques,  Feng  et  al. 
(1986a)  showed  the  relationships  for  the  first-  and  second-order  Lagran¬ 
gian  residual  currents,  which  has  already  been  discussed  in  Sections  2.2 
and  2.3.  Using  the  Lagrangian  residual  circulation  reported  by  Cheng 
(1983)  for  South  San  Francisco  Bay,  Feng  el  al.  (1986b)  compared  inter¬ 
preted  streamlines  with  observed  salinity.  According  to  their  theory 


for  the  zeroth  order  solution  for  salinity  transport,  the  isohaline  con¬ 


tours  of  salinity  should  be  identical  to  the  contours  of  the  streamlines 
for  the  steady-state  first-order  Lagrangian  residual  current.  The  ob¬ 
served  isohaline  contours  were  compared  qualitatively  with  the 
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streamline  contours.  The  agreement  was  relatively  good  considering  the 
computed  results  were  steady-state  and  did  not  include  stratification 
effects  (e.g.  their  model  was  2D  depth-averaged).  Feng  (1987)  made  a 
similar  comparison  of  salinity  contours  and  residual  currents  computed 
from  a  2D  depth-averaged  model  of  the  Bohai  Sea,  China.  Neither  Feng 
et  al.  (1986b)  nor  Feng  (1987)  attempted  to  use  their  residual  currents 
to  drive  a  salinity  (or  other)  transport  model. 

Orbi  and  Salomon  (1988)  also  used  a  2D  depth- averaged  model  to 
study  circulation  in  the  English  Channel.  Like  Cheng  and  Casulli  (1982) 
and  Cheng  (1983) ,  Orbi  and  Salomon  used  particle  tracking  throughout  the 
tidal  period  to  compute  Lagrangian  residuals.  They  also  found  a  spread 
of  the  residual  vectors,  or  tidal  current  phase  induced  mixing.  This 
system  is  relatively  nonlinear  ( k  -  0.25),  as  is  South  San  Francisco 
Bay.  Orbi  and  Salomon  also  did  not  use  residual  currents  to  drive  a 
transport  model. 

Najarian  et  al.  (1982  and  1984)  used  a  2D  laterally -averaged  model 
to  develop  intratidal  flow  fields  from  which  tidally- induced  residual 
currents  were  computed.  Eulerian  residuals  and  Stokes'  drift  were  used 
to  obtain  the  first-order  Lagrangian  residuals.  Najarian  et  al.  veri¬ 
fied  their  computed  residual  currents  by  comparing  model  results  with 
known  analytical  solutions  of  Lagrangian  residuals  in  a  homogeneous  es¬ 
tuary.  Model  experiments  were  performed  to  investigate  the  influences 
of  tidal  forcing,  density  gradients,  and  topographic  variations  on  mean 
Lagrangian  currents.  Najarian  et  al.  also  applied  the  model  t  O  thC 
Chesapeake  and  Delaware  Canal  to  assess  tide- induced  and  density- induced 
currents  in  the  canal.  However,  they  did  not  use  computed  residual  cur¬ 
rents  to  drive  a  transport  model. 
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Gomez-Reyes  (1989)  investigated  mechanisms  contributing  to  tidally 
driven  Lagrangian  residual  velocities.  He  expanded  to  second-order  the 
Euler-Lagrange  transformation  to  relate  Lagrangian  residual  currents  to 
the  spatially  and  time -varying  Euler ian  velocity  field.  The  Eulerian 
velocity  field  was  subjected  to  tidal  rectification  to  include  mean  vel¬ 
ocity,  M2,  and  M4  components.  The  second- order  approximations,  which 
are  referred  to  as  Lagrangian  drift  (Feng  et  al.  1986a),  include  terms 
that  are  functions  of  the  initial  release  times  of  particles.  Lagrangi¬ 
an  residuals  were  computed  from  the  results  of  a  depth -averaged,  numeri¬ 
cal  model  of  a  geometrically  simple,  shallow  bay  with  a  headland.  A 
grid  resolution  of  400  m  was  used,  and  k  was  on  the  order  of  0.10, 

These  results  were  compared  with  Lagrangian  residuals  computed  from  nu¬ 
merically  simulated  particle  trajectories.  Gomez-Reyes  concluded  that 
the  validity  of  the  Euler-Lagrange  transformation  depends  upon  the  de¬ 
gree  of  non-linearity  as  defined  by  the  ratio  of  the  local  tidal  excur¬ 
sion  to  the  local  length  scale  of  the  Eulerian  velocity  gradient.  For 
regions  in  the  interior  of  the  bay  where  the  non-linearity  ratio  was 
less  than  0.5,  the  first-order  approximation  (i.e.  Eulerian  residual 
plus  Stokes'  drift)  is  sufficient.  For  regions  close  to  the  headland 
where  the  non-linearity  ratio  was  between  0.5  and  1.0,  the  second-order 
approximation  is  necessary.  For  regions  immediately  next  to  the  head¬ 
land  where  the  ratio  is  greater  than  1.0,  Lagrangian  residual  currents 
should  be  computed  from  particle  trajectories.  Gomez-Reyes  did  not  use 
residual  currents  for  transport. 

The  literature  indicates  that  all  previous  applications  of  Lagran¬ 
gian  residual  currents  used  either  2D  depth-  or  2D  laterally- averaged 
hydrodynamic  model  results.  Also,  previous  residual  circulation  studies 
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typically  used  dynamic  steady-state  tidal  velocities  to  calculate  the 
Lagrangian  residual  circulation. 

None  of  the  calculated  Lagrangian  residual  circulation  fields  dis¬ 
cussed  above  were  used  to  drive  a  transport  (e.g.  water  quality)  model. 
As  mentioned  in  Section  2.1,  there  are  examples  of  the  use  of  Eulerian 
residuals  (arithmetic  averages  at  fixed  points  over  one  or  more  tidal 
cycles)  for  driving  salinity  and  water  quality  transport  models.  Also, 
intratidal  velocities  have  been  used  for  intratidal  transport,  as  in  the 
salinity  simulation  by  Cheng  et  al.  (1984)  and  the  water  quality  simula¬ 
tions  by  Hall  (1989) .  There  are  not  any  known  examples  of  the  use  of 
Lagrangian  residual  circulation  for  driving  a  transport  model.  Also, 
there  are  not  any  known  examples  of  computing  time-varying,  3D,  Lagran¬ 
gian  residual  circulation  from  an  intratidal  hydrodynamic  model. 

2 . 5  CHAPTER  SUMMARY 

Lagrangian  residual  currents  should  be  used  to  advect  water  masses 
for  long-term,  intertidal  transport  models.  The  Lagrangian  residual 
transport  equation,  which  uses  Eulerian  field  variables,  has  the  form  of 
Equation  2.10.  First-order  Lagrangian  residual  currents,  which  are  the 
sum  of  the  Eulerian  residual  velocity  and  the  Stokes'  drift,  should  be 
adequate  for  weakly  nonlinear  systems,  such  as  Chesapeake  Bay.  Formula¬ 
tions  for  Stokes'  drift  shown  in  Equations  2.17  and  2.18  can  be  used  to 
ensure  mass  conservation  for  proper  use  in  a  transport  model.  A  litera¬ 
ture  search  revealed  that  there  are  no  previous  studies  of  using  Lagran¬ 
gian  residual  circulation  to  drive  an  intertidal  transport  mo^el. 


CHAPTER  3 


COMPUTATIONAL  METHODOLOGY 

Descriptions  of  the  HM,  WQM,  and  an  interface  processor  within  the 
HM  are  presented  in  this  chapter.  The  interface  processor  computes  and 
outputs  residual  currents  and  other  HM  information  for  input  to  the  WQM. 
The  mathematical  adaptation  and  numerical  implementation  of  the  residual 
computations  within  the  interface  processor  are  also  discussed. 

3.1  MODELING  FRAMEWORK  DESCRIPTION 

This  section  of  the  chapter  describes  the  three  components  of  the 
modeling  framework,  the  HM,  WQM,  and  interface  processor.  The  basic 
features  of  the  interface  processor  are  introduced  within  this  section, 
while  the  details  of  the  processor  are  presented  in  the  following  sec¬ 
tions  (i.e.  Sections  3. 2 -3. 4). 

3.1.1  Hvdrodvnamic  Model 

The  hydrodynamic  model  used  in  this  study  is  referred  to  as  Cur¬ 
vilinear  Hydrodynamics  in  Three  Dimensions  (CH3D) .  CH3D  is  a  finite 
difference  model  for  calculating  free  surface,  time -varying,  three- 
dimensional  currents,  temperature,  and  salinity  in  surface  waters  (e.g, 
estuaries,  lakes,  and  coastal  embayments) .  The  model  originated  with 
Cartesian  horizontal  coordinates  and  a  stretched  vertical  coordinate 
(Sheng  1983  and  1984) ,  but  was  subsequently  modified  for  curvilinear 
(boundary- fitted)  coordinates  in  the  horizontal  plane  (Sheng  1986a  and 
1986b) .  The  model  was  modified  for  use  on  Chesapeake  Bay  by  Johnson  et 
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al.  (1989)  to  allow  a  fixed  Cartesian  vertical  coordinate.  Therefore, 
the  version  used  in  this  study  has  boundary- fitted  coordinates  in  the 
horizontal  plane  and  fixed  layers  in  the  vertical  direction. 

Boundary- fitted  coordinates  provide  a  generalized  method  of  mapping 
irregular  geometry  without  requiring  orthogonality  (Thompson  and  Johnson 
1985) .  The  physical  and  transformed  grids  used  for  the  Chesapeake  Bay 
model  are  shown  in  Figure  3.1.  The  distances  between  grid  lines  in  the 
transformed  grid  are  assigned  values  of  unity  for  convenience  in  the 
solution  algorithms. 

CH3D  solves  the  3D  equations  of  continuity,  momentum,  and  conserva¬ 
tion  of  salt  and  temperature.  Salt  and  temperature  are  related  to  den¬ 
sity  through  an  equation  of  state.  The  model  can  also  transport  a  con¬ 
servative  tracer.  The  hydrostatic  pressure  assumption  is  used  for  the 
vertical  momentum  equation.  The  baroclinic  pressure  effects  are  re¬ 
tained  in  the  momentum  equations  through  the  coupling  of  salt  and/or 
temperature  transport  with  momentum.  All  physical  processes  impacting 
estuarine  circulation  are  included  in  the  model,  such  as  tides,  wind, 
density  effects,  freshwater  inflows,  the  earth's  rotation,  and  tur¬ 
bulence.  Several  higher  order  closure  schemes  are  available  for  solving 
the  vertical  eddy  viscosity  and  diffusivity.  A  simplified  second-order 
vertical  turbulence  model,  which  is  based  on  the  assumption  of  local 
equilibrium  of  turbulence,  was  selected  for  use  on  Chesapeake  Bay 
(Johnson  et  al.  1989). 

The  governing  equations  are  transformed  into  boundary- fitted  co¬ 
ordinates  (^ , »7 ) .  Both  the  planform  Cartesian  coordinates  (x,y)  and 
velocities  (u,v)  are  transformed  into  the  curvilinear  system  such  that 
the  velocity  components  are  normal  to  the  (£,>?)  coordinate  lines.  This 
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Figure  3.1(a).  Physical  grid  of  Chesapeake  Bay 
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Figure  3.1(b).  Transformed  grid  of  Chesapeake  Bay 
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is  accomplished  by  defining  the  Cartesian  velocities  in  terms  of  contra- 
variant  components  (U,V).  The  governing  equations  are  also  nondimen- 
sionalized  within  CH3D.  Basic  concepts  concerning  the  transformations 
and  the  nondimens ional  relationships  are  presented  in  Appendix  A.  The 
details  of  CH3D's  governing  equations  can  be  found  in  Sheng  (1986a  and 
1986b)  for  stretched  vertical  coordinates  and  Johnson  et  al.  (1989)  for 
fixed  vertical  layers. 

The  transformed  equations  are  solved  on  the  transformed  grid  for 
the  dependent  variables  U,  V,  W  (vertical  physical  velocity),  S  (sali¬ 
nity),  T  (temperature),  and  f  (water  surface  displacement).  Although 
the  dependent  variables  are  nondimens ional  quantities,  the  asterisks  are 
omitted  here  for  convenience. 

The  solutions  for  the  dependent  variables  are  accomplished  through 
external  and  internal  modes.  The  external  mode  solves  for  f  and  ver¬ 
tically  integrated  contravariant  flow  per  unit  width  from  the  trans¬ 
formed,  vertically  integrated,  horizontal  momentum  and  continuity  equa¬ 
tions.  All  terms  in  the  continuity  equation  and  the  water  surface  slope 
terms  of  the  momentum  equations  are  treated  implicitly  and  weighted  be¬ 
tween  new  and  old  time-levels.  The  resulting  equations  are  factored 
such  that  a  £-sweep  followed  by  an  q-sweep  yields  the  solution  at  the 
new  time -level. 

The  velocities  U  and  V  are  solved  from  the  transformed,  layer- 
averaged,  horizontal  momentum  equations;  W  is  solved  from  the  trans¬ 
formed,  layer -averaged  continuity  equation;  and  S  and  T  are  solved  from 
the  transformed,  layer -averaged  mass  and  heat  conservation  equations. 

The  solutions  of  these  equations  constitute  the  internal  mode.  The 
vertical  diffusion  terms  of  all  the  internal  mode  equations  and  the 


34 


bottom  friction  and  water  surface  slope  terms  of  the  momentum  equations 
are  treated  implicitly.  Water  surface  deviations  from  the  external  mode 
are  used  in  the  water  surface  slope  terms  of  the  internal  mode.  After  U 
and  V  are  computed,  they  are  adjusted  to  ensure  mass  conservation  by 
forcing  the  sum  of  the  velocities  over  the  vertical  to  equal  the  verti¬ 
cally  integrated  unit  flow.  Finally,  contravariant,  nondimens ional 
quantities  are  converted  to  physical,  dimensional  quantities.  Veloci¬ 
ties  in  Cartesian  coordinates  are  also  computed  for  plotting  velocity 
vectors  on  the  physical  grid. 

3.1.2  Water  Quality  Transport  Model 

The  water  quality  model  used  for  the  Chesapeake  Bay  study  and  this 
study  was  recently  developed  at  the  US  Army  Engineer  Waterways  Experi¬ 
ment  Station  (WES).  The  model  is  based  on  the  integrated  compartment 
method,  ICh  (i.e.  the  mass  transport  equation  is  applied  in  integrated 
form  to  control  volumes,  or  boxes).  The  ICM,  or  box  model  approach,  was 
followed  from  the  USEPA  Water  Quality  Analysis  Simulation  Program  (WASP) 
(DiToro  et  al.  1983  and  Ambrose  et  al.  1986)  to  facilitate  coupling  with 
various  hydrodynamic  models.  Also,  the  ICM  facilitates  overlaying  the 
HM  grid  with  a  coarser  WQM  grid  (Bird  and  Hall  1988  and  Hall  et  al. 

1988) .  The  ICM  has  been  used  by  others  for  different  numerical  modeling 
applications,  such  as  solution  of  the  Navier-Stokes  equations  by  Yeh 
(1981). 

The  basis  for  the  ICM  is  conversion  of  the  mass  conservation  equa¬ 
tion  for  an  infinitesimal  point  (Equations  2.1)  into  a  finite  control 
volume  (i.e.  integrated  compartment)  form  from  which  it  was  originally 
derived.  Using  the  integral  theorems  of  vectors  and  replacing  integrals 
over  surfaces  with  summations  over  compartment  (segment)  faces,  Equation 
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2.1  can  be  cast  (Yeh  1981)  as 

3(C,V.)  „  „  D,,An 

"  £  QijCU  +  l  (Gj  -  ci >  ±  X  (SOURCES/SINKS).  (3.1) 

oc  j  J  J  O  O 

where 

i  -  segment  index 

j  -  segment  index  of  adjoining  segment 
Vi  -  segment  i  volume 
C£  -  segment  i  concentration 

Qij  -  flow  to  (positive)  or  from  (negative)  segment  i 
from/to  segment  j 

Cij  -  concentration  at  interface  of  segment  i  and  j 

Djj  -  eddy  diffusion  coefficient  for  the  ij  interface 

Aij  -  facial  area  of  the  ij  interface 

Lij  -  mixing  length  (box  lengths)  between  segments  i 
and  j 

(SOURCES/S INKS )L  -  rate  of  change  of  mass  in  segment  i  from 

m  various  sources  and/or  sinks  m  (loads,  kinetic 
transformations ,  etc . ) 


It  should  be  noted  that  flows ,  rather  than  velocities ,  are  used 
with  the  ICM.  In  the  WASP  code,  the  finite  difference  approximations 
for  the  tiarms  of  Equation  3.1  are  the  same  for  all  directions.  The  WASP 
code  uses  central  differencing  for  diffusion.  The  user  can  select 
either  of  two  interpolation  methods  for  C^j ,  which  results  in  either 


backward  (upwind)  or  central  differencing  for  advection.  With  the  for¬ 


ward  in  time  (explicit)  solution  scheme  of  WASP,  central  differencing 
can  often  result  in  numerical  instabilities  (Roache  1972  and  Ambrose  et 
al.  1986).  Use  of  the  upwind  differencing  for  advection  can  result  in 
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excessive  numerical  diffusion  (Roache  1972  and  Ambrose  et  al.  1986). 
Explicit  solution  schemes  for  the  vertical  dimension  can  result  in  over¬ 
ly  restrictive  time  steps  when  simulating  periods  of  high  vertical  mix¬ 
ing  (Sheng  1983) . 

The  structure  of  the  WQM  retained  the  ICM,  but  the  solution  schemes 
were  rebuilt  to  address  the  concerns  expressed  in  the  previous  para¬ 
graph.  The  WQM  distinguishes  between  the  horizontal  and  vertical  direc¬ 
tions.  Horizontal  advective  fluxes  are  normally  several  orders  of  mag¬ 
nitude  greater  than  diffusive  fluxes  in  surface  waters;  thus,  it  is 
desirable  to  accurately  resolve  advective  fluxes  without  introducing  nu¬ 
merical  diffusion  that  is  greater  than  the  physical  diffusion.  The  use 
of  higher-order  advection  schemes  have  dramatically  reduced  or  eliminat¬ 
ed  the  concerns  associated  with  numerical  diffusion/dissipation  in 
Eulerian  transport  models.  The  Quadratic  Upstream  Interpolation  for 
Convective  Kinematics  with  Estimated  Streaming  Terms  (QUICKEST)  scheme 
(Leonard,  1979),  which  is  explicit,  upstream  weighted  and  third-order 
accurate  in  space,  was  implemented  (Chapman  1988)  for  horizontal  advec¬ 
tion  in  the  WQM.  Either  QUICKEST  or  upwind  differencing  are  user  op¬ 
tions.  CH3D  also  has  the  option  to  use  QUICKEST  or  upwind  for  horizon¬ 
tal  advection  of  mass  (e.g.  salinity). 

The  solutions  for  the  horizontal  and  vertical  transport  use  a  split 
scheme.  First,  horizontal  advection  and  diffusion  for  all  cells  are 
computed  explicitly  with  a  two  time  level  forward  time  difference  to 
provide  a  provisional  update  for  constituent  concentrations.  The  trans¬ 
port  solution  for  constituent  concentrations  at  the  next  time  level  is 
completed  by  updating  the  provisional  concentrations  with  vertical  ad¬ 
vection  and  diffusion.  The  Crank-Nicolson  scheme,  which  is 
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unconditionally  stable  (Roache  1972) ,  is  used  for  vertical  advection  in 
the  WQM.  The  user  can  select  the  time  weighting  factor  between  0.5  and 
1.0,  where  a  value  of  1.0  results  in  a  fully  implicit,  central  dif¬ 
ference  scheme.  A  fully  implicit,  central  difference  scheme  is  used  for 
vertical  diffusion.  The  vertical  sweeps  are  conducted  column  by  column 
with  a  tridiagonal  matrix  solved  for  cell  concentrations  within  each 
column.  This  approach  removes  the  stability  restriction  associated  with 
small  vertical  box  lengths. 

The  transport  solution  procedures  are  conducted  for  each  water 
quality  constituent.  If  a  water  quality  constituent  is  reactive  (i.e. 
non-conservative) ,  then  the  concentration  changes  resulting  from  the 
various  sources/sinks  are  added. 

The  WQM  allows  time -varying  boundary  conditions  and  hydrodynamic 
updates.  Also,  the  user  can  specify  a  constant  model  time  step  or  se¬ 
lect  the  autostepping  feature,  which  automatically  adjusts  the  time  step 
to  satisfy  the  horizontal  flow  stability  restriction.  This  feature  was 
included  to  take  advantage  of  potentially  larger  time  steps  during  low 
flow  periods  of  the  simulation. 

3.1.3  Interface  Processor 

The  interface  processor  couples  the  HM  and  WQM  computational  grids 
and  processes  hydrodynamic  information  into  WQM  input  data.  The  inter¬ 
face  processor  was  developed  as  subroutines  within  the  HM.  Therefore, 
the  hydrodynamic  information  for  the  WQM  is  processed  and  stored  while 
the  HM  is  executing. 

Coupling  of  the  HM  and  WQM  grids  requires  generation  of  map  files, 
which  set  up  a  correspondence  between  the  HM  and  WQM  grid  formats. 
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Additionally,  time -invariant  HM  geometric  information  is  required  to 
compute  box  lengths,  voltues,  and  facial  areas  for  use  by  the  WQM. 

Processing  of  the  time-varying  hydrodynamic  information  into  WQM 
input  data  can  be  accomplished  in  either  of  two  modes,  intratidal  or 
intertidal.  The  intratidal  mode  involves  processing  the  hydrodynamics, 
produced  at  intervals  on  the  order  of  minutes  (e.g.  five  minutes),  into 
WQM  input  at  about  one-  or  two-hour  intervals.  The  intratidal  mode  sim¬ 
ply  requires  temporal  averages  of  the  hydrodynamics  (flows  and  vertical 
diffusivities) .  The  intertidal  mode  involves  processing  hydrodynamics 
into  WQM  input  at  tidal-period  intervals  or  greater,  thus,  reducing  WQM 
input  data  storage  requirements  by  an  order  of  magnitude.  Intertidal 
processing  requires  computation  of  the  Eulerian  residuals  and  Stokes' 
drifts.  For  both  modes,  processed  flows  and  diffusivities  are  output  in 
a  format  compatible  with  the  WQM  following  appropriate  scaling.  Scaling 
accounts  for  the  fact  that  the  contravariant  velocities  in  the  HM  are 
both  nondimens ional  and  defined  on  a  transformed  boundary- fitted  grid. 

The  WQM  can  accept  either  intratidal  or  intertidal  input.  The  only 
difference  is  that,  for  the  intertidal  mode,  Eulerian  residuals  and 
Stokes'  drifts  that  are  input  to  the  WQM  are  added  together  to  produce 
the  Lagrangian  residuals,  which  are  used  in  the  WQM  transport  equation. 
Only  Eulerian  residuals  are  used  for  the  intratidal  mode. 

3.2  COMPUTATIONAL  PROCEDURES  OF  THE  INTERFACE  PROCESSOR 

The  interface  processor  (i.e.  subroutines  within  the  HM)  performs 
three  basic  tasks:  1)  maps  the  WQM  grid  to  the  HM  grid;  2)  outputs 
time -invariant  HM  geometric  information  required  by  the  WQM;  and  3)  out¬ 
puts  time -varying  HM  transport  information  used  for  WQM  transport 
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computations.  The  procedures  and/or  computations  for  each  process  are 
described  below. 

3.2.1  Grid  Mapping 

The  WQM  grid  must  be  mapped  to  the  HM  grid  because  the  two  models 
have  distinctly  different  numerical  frameworks.  The  HM  uses  an  ijk 
coordinate  system  that  corresponds  to  the  transformed  plane  on  which  the 
equations  of  motion  are  solved.  The  WQM  uses  a  sequential  cell  number¬ 
ing  configuration.  Numbering  progresses  along  the  horizontal  plane 
starting  with  the  surface  layer.  Each  surface  cell  can  have  multiple 
cells  (layers)  beneath  it.  Therefore,  HM  cells  are  referred  to  by  i,j,k 
indices,  whereas  WQM  cells  are  referred  by  cell  (box)  numbers.  An  ex¬ 
ample  of  HM  and  WQM  cell  numbering  for  a  simple  3x3x3  grid  is  shown  in 
Figure  3.2. 

Although  the  box  numbering  concept  is  cumbersome,  it  does  allow 
flexibility  in  configuring  the  WQM  grid.  For  example,  it  is  possible  to 
overlay  the  HM  grid  with  a  coarser  WQM  grid,  as  was  done  in  a  study  of 
Los  Angeles  -  Long  Beach  Harbors  by  Hall  (1989) .  Grid  overlaying  was 
not  used  in  this  research;  the  HM  and  WQM  use  the  same  grid  configura¬ 
tion  and  density  for  all  developments  and  tests. 

The  mapping  procedure  begins  by  creating  two  map  files  that  are 
read  into  the  interface  subroutine,  WQMOUT.  One  of  the  map  files, 
FILE94,  declares  time -averaging  parameters  in  addition  to  mapping  infor¬ 
mation.  FILE94  specifies  the  number  of  boxes  in  the  surface  layer 
(NSB) ,  the  averaging  interval  of  the  HM  time  iterations  (NAVG) ,  and  the 
HM  iteration  number  on  which  time  averaging  is  to  start  (ITWQS) .  The 
variables  in  parentheses  are  FORTRAN  variables  used  in  the  processor, 
which  is  presented  in  Appendix  B.  FILE94  also  specifies  the  cell  or  box 
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gure  3.2.  Cell  numbering  schemes  for  the  hydrodynamic 
and  water  quality  models 
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number  (NB)  for  each  WQM  surface  cell  and  the  starting  and  ending  HM 
horizontal  plane  grid  lines  (I,J)  that  define  that  surface  box 
(IFIRST(SB),  ILAST(SB),  JFIRST(SB) ,  AND  JLAST(SB)).  For  a  one-to-one 
grid  correspondence,  the  starting  and  ending  grid  line  indices  differ  by 
one.  These  grid  line  specifications  are  only  necessary  for  the  surface 
layer,  since  all  other  layers  fall  within  the  same  horizontal  grid  as 
the  surface  layer. 

Recall  from  Section  3.1.2  that  transport  in  the  WQM  is  handled  on  a 
cell  by  cell  basis  where  flows  and  dif fusivities  are  specified  on  cell 
faces.  Additionally,  the  WQM  does  not  have  separate  directional  indices 
for  the  horizontal  plane.  Therefore,  sufficient  information  must  be 
provided  to  identify  horizontal  linkage  of  flow  faces  and  boxes. 

FILE95,  which  is  used  by  both  WQMOUT  and  the  WQM,  provides  mapping  in¬ 
formation  for  relating  horizontal  flow  faces  and  WQM  boxes.  FILE95  spe¬ 
cifies  the  number  of  horizontal  flow  faces  for  the  surface  layer  (NHQF) 
and  the  total  number  of  horizontal  flow  faces  for  all  layers  (NHQFT) . 
FILE95  also  specifies  for  each  horizontal  flow  face  (F-l, NHQFT)  the  fol¬ 
lowing  : 

FN  -  horizontal  flow  face  number 

QD(F)  -  flow  direction  code  (QD=1  for  U,  2  for  V  directions) 

IL(F)  -  WQM  cell  to  the  left  (upstream  direction)  of  cell  IQ 
(see  Figure  3.3) 

IQ(F)  -  WQM  cell  to  the  left  (upstream  direction)  of  flow  face  F 

JQ(F)  -  WQM  cell  to  the  right  (downstream  direction)  of  flow  face  F 

JR(F)  -  WQM  cell  to  the  right  (downstream  direction)  of  cell  JQ 

KP(F)  -  index  of  HM  grid  line  corresponding  to  flow  face  F 


QD(F)=2 
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Figure  3.3.  Nomenclature  for  map  file  information 
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KF(F) ,  KL(F)  -  range  of  HM  horizontal  cells  for  spatial  averaging 
(overlay)  of  flows  for  WQM.  For  1:1  overlay,  KF-KL.  When 
KP  -  I  index  (QD-1) ,  KF,  KL  -  J  index.  When  KP  -  J  index 
(QD-2) ,  KF,  KL  -  I  index. 

KZ(F)  -  K  index  (layer)  of  HM  for  horizontal  flow  face  F. 

NHQF,  NHQFT,  QD,  KP,  KF,  and  KZ  are  used  within  WQMOUT  and  FN,  XL,  IQ, 

JQ,  JR,  and  KZ  are  used  within  the  WQM  to  relate  HM  flows  to  the  WQM. 

Mapping  flows  in  the  vertical  direction  is  much  more  straight¬ 
forward  because  the  WQM  has  sense  of  the  vertical  direction.  The  HM 
contains  the  array  KM(I,J)  which  defines  the  bottom  layer  of  every  IJ 
column  as  shown  in  Figure  3.4.  The  vertical  flows  and  diffusivities 
computed  by  the  HM  are  converted  and  stored  by  WQMOUT  as  2D  arrays 
(SB,K) ,  where  SB  is  the  surface  box  number  and  K  is  the  layer  between  KM. 
and  KMAX.  Therefore,  mapping  of  vertical  flows  and  diffusivities  is 
accomplished  by  looping  over  the  layers  for  each  surface  cell  while 
sweeping  all  the  surface  cells . 

3.2.2  Geometric  Information 

Time -invariant  geometric  information  is  provided  to  the  WQM  by  the 
interface  subroutine  WQMOUT.  This  information  consists  of  the  follow¬ 
ing  with  processor  variables  noted  in  parentheses: 

AZ  (DELTAZ)  -  layer  thicknesses  for  layers  beneath  the  surface 

layer,  m.  All  are  the  same  thickness  and  constant  with  time. 

AZra  (DELTAZM)  -  nominal  layer  thickness  for  the  surface  layer,  m. 
Total  surface  layer  thickness  varies  with  time  and  includes 
the  deviation  of  the  water  surface  from  DELTAZM. 
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the  hydrodynamic 
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Ag  (BSAREA(SB))  -  plan  surface  area  of  each  column  (surface  box), 

m2. 

(BL(SB,1)),  (BL(SB,2))  -  horizontal  cell  dimensions  (box 

lengths)  in  the  £  and  t\  directions,  respectively,  m. 

AH^  and  AH^  (FAREA(F))  -  cell  facial  area  for  horizontal  diffusivi- 
ty  in  the  £  and  17  directions,  respectively,  m2.  FAREA  is  con¬ 
stant  with  time  except  for  the  surface  layers;  thus,  FAREA  for 
the  surface  layer  is  computed  and  output  when  time  averaging 
starts  (ITWQS) . 

Vc  (CVOL(SB))  -  cell  volumes  for  all  cells  beneath  the  surface  lay¬ 
er,  m^.  All  are  time -invariant  and  the  same  volume  within  a 
column  since  layer  thicknesses  beneath  the  surface  layer  are 
uniform  thickness  and  constant  with  time. 

Vs  (CVOLS(SB))  -  cell  volumes  for  all  cells  in  the  surface  layer, 

m2 .  CVOLS  vary  with  time ;  thus ,  CVOLS  are  computed  and  output 
when  time  averaging  starts  (ITWQS) . 

The  layer  thicknesses  and  box  lengths  are  used  in  the  advection  and 
diffusion  calculations  of  the  WQM.  Cell  face  areas  (BSAREA  and  FAREA) 
are  used  only  in  the  diffusion  calculations.  Cell  face  flows  used  in 
advection  calculations  are  computed  within  WQMOUT;  thus,  facial  areas 
are  not  needed  for  advection.  Initial  ceil  volumes  are  required  for 
initial  conditions  in  the  WQM. 

The  dimensional  scaling  and  conversion  from  contravariant  to  physi¬ 
cal  components  (taking  coordinate  transformations  into  account)  for  the 
geometric  quantities  are  shown  below  (also  refer  to  Appendix  A) . 
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AZ  and  AZm  -  Since  the  vertical  dimension  is  not  transformed  and 
layer  thicknesses  exist  in  dimensional  and  nondimensional  form,  AZ 
and  AZm  need  only  to  be  converted  from  centimeters  to  meters. 


A.  - 

(3.2) 

where  g01/2  is  the  square  root  of  the  determinant  of  the  metric  ten¬ 
sor  (i.e.  Jacobian  of  the  transformation),  located  at  the  cell  cen¬ 
ter.  Xr  is  a  dimensional  scale  factor  for  horizontal  lengths. 

(3.3) 

“  4®22  Xr 

(3.4) 

where  gu  and  g22  are  the  metric  coefficients  which  scale 

transformation  in  the  £  and  rj  directions,  respectively. 

gu  and  g22  are  taken  at  the  cell  center. 

the  grid 

Values  for 

AH*  “  AZ*  Zr  ^  X, 

(3.5) 

AH„  -  AZ*  Zr  ^  Xr 

(3.6) 

where  AZ  is  the  nondimensional  layer  thickness.  The  metric  coef¬ 
ficients  are  defined  at  the  respective  U  and  V  faces.  Zr  is  a 
dimensional  scale  factor  for  vertical  lengths. 

Vc  -  AZ  Ae 


(3.7) 


47 


V8  -  (A2^  +  r*  rr)  a,  (3.8) 

where  f  is  the  nondimens ional  distance  from  the  top  of  AZg,  to  the 
water  surface.  is  the  dimensional  scale  factor  for  the  water 
surface  deviation. 

All  HM  dimensional  length  scales  are  in  centimeters  and  are  converted  to 
meters  for  the  WQM  prior  to  output. 

3.2.3  Transport  Information 

Transport  information,  which  consists  of  flows  for  all  cell  faces 
and  vertical  eddy  diffusivity  coefficients,  is  processed  within  the  HM 
and  output  for  use  in  the  WQM.  Transport  information  is  averaged 
throughout  each  HM  averaging  interval  (NAVG)  and  output  to  a  disk  file 
at  the  end  of  the  interval.  Additionally,  HM  elapsed  simulation  time, 
surface  cell  volumes  (CVOLS) ,  and  surface  cell  facial  areas  for  horizon¬ 
tal  diffusive  fluxes  (FAREA)  are  computed  and  output  at  the  end  of  each 
time -averaging  interval.  Surface  cell  volumes  are  required  in  the  WQM 
to  compute  water  surface  elevation  and  to  check  volume  and  mass  bal¬ 
ances  . 

Horizontal  flow  in  the  £  direction,  Q^  (HQ(F)),  is  a  dimensiona1 
physical  component  computed  from 

Qf-U*Ur  ^XtAZ*Zr  (3.9) 

where 

U*  -  dimensionless  contravariant  velocity  in  the  £  direction 
^  “  Jacobian  defined  at  the  U  flow  face 


Conversion  from  contravariant  to  physical  components  results  in  flows 
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that  are  still  normal  to  transverse  grid  lines  (i.e.  cell  faces)  as 
required  by  the  WQM. 

Horizontal  flow  in  the  rj  direction  is  computed  by 

%  "  Ur  ^  X*  AZ*  Zr  (3.10) 


where 

V*  -  dimensionless  contravariant  velocity  in  the  r)  direction 


The  Jacobian  is  defined  at  the  V  flow  face. 

Vertical  flow,  Qt  (ZQ(SB,K)),  is  computed  by 

Q*“W*Ur  ^XrZr  0.11) 


where 

W*  -  dimensionless  physical  velocity  in  the  vertical  direction 
The  Jacobian  is  defined  at  the  cell  center. 

Vertical  diffusivity,  Dt  (DIFZ(SB,K)) ,  is  computed  according  to 

o.  -  \r  (3.12) 


where 

Gb  -  nondimens ional  vertical  eddy  diffusivity  computed  by  the  HM 
-  scale  factor  for  dimensional  vertical  eddy  diffusivity 

It  is  not  necessary  to  compute  and  output  horizontal  eddy  diffusivities 
since  the  HM  uses  a  constant  value  for  this  parameter,  as  does  the  WQM. 
Horizontal  eddy  diffusivity  is  relatively  unimportant  in  advection 
dominated  systems  such  as  estuaries  (Feng  et  al.  1986b  and  Bedford 
1985),  with  possible  exceptions  very  near  boundaries. 
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The  arithmetic  time  average  of  a  variable  can  be  computed  from 


<t> 


(3.13) 


where 

N  -  number  of  intervals  averaged 
Equation  3.13  is  implemented  in  WQMOUT  by  an  equivalent  form, 

\  “  ?n-l  +  ¥  (3>14) 

Therefore,  each  transport  variable  is  updated  each  HM  time  step.  When  a 
counter  in  WQMOUT  reaches  N  (NAVG) ,  the  end  of  the  averaging  interval 
has  been  attained,  and  the  time -averaged  variables  are  output. 

Time-averaged  flows  consist  of  Eulerian  residuals  and  Stokes' 
drifts.  Equations  3.9-3.11  are  used  along  with  Equation  3.14  to  compute 
the  Eulerian  residuals.  Equations  3,12  and  3.14  are  used  to  compute  the 
time  average  of  the  vertical  diffusivities.  The  computations  for  the 
Stokes'  drift  flows  are  presented  in  the  next  section. 

3.3  IMPLEMENTATION  OF  STOKES'  DRIFT  COMPUTATIONS 

Equations  2.17  and  2.18  are  the  basic  equations  used  for  computing 
the  Stokes'  drift.  The  approach  used  is  to  numerically  implement  Equa¬ 
tions  2.17  and  2.18  within  a  subroutine  in  the  HM  using  the  nondimen- 
sional,  contravariant  velocities  as  they  are  computed  by  the  HM  through¬ 
out  the  averaging  interval.  At  the  end  of  the  averaging  interval,  the 
computation  for  Stokes'  drift  is  completed,  the  velocities  are  converted 
to  dimensional,  physical  Stokes'  flows,  and  the  time-averaged  Stokes' 
flows  are  output. 
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Within  an  averaging  interval,  the  mean  (tidally  averaged)  veloci¬ 
ties  and  the  deviations  about  the  mean  velocities  (i.e.  U' ,  V' ,  W' , 
where  nondimens ional  variables  are  understood  and  the  asterisks  have 
been  dropped  for  convenience)  are  not  known.  It  appears  that  Equation 
2.18  can  not  be  implemented  while  the  HM  is  executing  without  storing 
the  entire  velocity  field  for  every  time  step  throughout  the  averaging 
interval.  However,  the  velocity  deviation  relationships  of  Equation 
2.18  can  be  obtained  from  mean  quantities.  Consider  the  Bz  term  of 
Equation  2.18  as  an  example.  This  term  can  be  written  in  a  shortened 
notation  as 

Bz  -  U'  /  V'  dt  -  U'  tj'  (3.15) 

where 

q'  -  cumulative  displacement  vector  resulting  from  the  velocity 
deviation  V' 

Substituting  the  decomposition  relationship  of  Equation  2.2,  expanding 
Equation  3.15,  taking  averages,  and  cancelling  terms  results  in 

\  -  \Tri  -  U~t  V  (3.16) 

where  t  is  the  elapsed  time  in  the  averaging  interval.  Therefore,  it  is 
possible  to  compute  the  B  terms  (vector  potentials) ,  with  minimal  stor¬ 
age  requirements,  by  accumulating  mean  quantities  during  the  averaging 
period  (using  Equation  3.14)  and  executing  equation  3.16  (or  similar 
equations  for  the  other  B  terms)  at  the  end  of  the  averaging  period. 

The  expanded  form  of  Equation  2.17  is 

83z  3By 
Us  “  ~dy  ■  lz 


(3.17) 
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V  “s 

«  dz  dx 


(3.18) 


w 

«  3x  3y 


(3.19) 


where  all  velocities  are  nondimensional  and  contravariant.  Although  the 
coordinates  in  Equations  3.17-3.19  should  be  transformed  coordinates 
($,rj),  Cartesian  variables  (x,y)  have  been  temporarily  interchanged  to 
avoid  confusion  with  the  displacement  variables  (£,ij,6)  discussed  below. 
By  locating  the  above  vector  potentials  as  shown  in  Figure  3.5,  the  spa¬ 
tial  derivatives  for  each  Stokes'  velocity  in  Equations  3.17-3.19  con¬ 
veniently  provide  second-order  central  differences. 

Spatial  averaging  of  velocities  and  their  integrated  displacements 
is  required  to  compute  the  vector  potentials  for  the  locations  shown  in 
Figure  3.5.  The  integrated  displacements  are 


(  -  /  U  dt 

(3.20) 

T)  -  f  v  dt 

(3.21) 

S  -  J  W  dt 

(3.22) 

where  the  above  velocities  are  defined  as  shown  in  Figure  3.5.  Two 
point  averages  between  adjacent  cells  are  used  to  obtain  the  B  terms  as 
follows : 


( 


Vi,j,k  +  Vl,J,k+l 


.k 


) 


(3.23) 
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n  fWi,J,k  +  Wi-l.J.O  +  ^i.J.k+ll 

\j,k  l - 2 - J  l - 2 - J  <3-24> 

\  j  k  “  (UiJ,1C  +2Ul,J~1,k)  (■—'*  +2l?i‘1,J,k)  (3.25) 


The  Bx  and  By  terms  at  the  water  surface  are  set  to  zero  to  keep 
the  formulation  conservative.  There  can  not  be  a  Stokes'  drift,  Ws,  and 
mass  flux  through  the  water  surface  in  a  mass  conserving  transport  mod¬ 
el.  The  vector  potentials  (i.e.  B  terms)  are  automatically  computed  as 
zero  along  solid  boundaries  with  the  above  formulations.  However,  at 
corner  points  of  solid  boundaries  that  protrude  into  the  flow  field,  it 
is  possible  to  compute  a  vector  potential  that  should  not  exist.  Such 
corners  can  exist  in  plan  and  elevation  views,  and  the  B  terms  must  be 
set  to  zero.  Likewise,  at  river  inflow  boundaries,  the  vector  poten¬ 
tials  are  set  to  zero.  At  tidal  boundaries,  the  B  terms  can  be  computed 
one  row  or  column  inside  the  grid  from  the  HM  tidal  boundary  (i.e.  where 
water  surface  is  specified).  Eulerian  and  Stokes'  flows  can  be  computed 
along  the  inner  grid  line.  Therefore,  along  the  tidal  boundary,  the  WQM 
grid  should  start  one  row  or  column  within  the  HM  grid. 

Equations  3.17-3.19  use  the  vector  potentials  to  compute  the 
Stokes'  velocities.  However,  it  is  Stokes'  flows  (volume/time) ,  rather 
than  velocities,  that  are  required  for  the  WQM.  To  keep  the  Stokes' 
flows  in  the  mass  conserving  form  of  Equation  2.17,  the  vector  poten¬ 
tials  must  be  converted  from  a  velocity  related  quantity  to  a  flow 
related  quantity.  This  conversion  is  accomplished  in  boundary- fitted 
coordinates  by  multiplying  the  vector  potentials  (i.e.  B  terms,  such  as 
Equation  3.16)  by  the  Jacobian  of  the  transformation, 
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Ai  -  4^  B* 


(3.26) 


where  the  subscript  i  represents  the  coordinate  directions.  Spatial 
averaging  of  the  Jacobian  is  required  for  the  Bz  terras  since  the  Jaco¬ 
bian  is  not  explicitly  defined  in  CH3D  at  intersections  of  grid  lines. 

The  A  terms  from  Equation  3.26  are  used  in  Equations  3.17-3.19 
rather  than  the  B  terms,  and  the  results  are  scaled  to  produce  formula 
tions  for  the  physical,  dimensional  Stokes'  flows  in  the  £,ri,z  direc¬ 
tions  (reverting  back  to  true  transformed  coordinate  variables).  The 
physical,  dimensional  Stokes'  flows  are,  thus,  computed  from 


AZ*  Ur  Xr  Zr  R0 
AZ*  Ur  Xr  Zr  R0 
Ur  Xr  Zr  R„ 


(3.27) 

(3.28) 

(3.29) 


where  RQ  is  the  Rossby  number.  The  Rossby  number  is  defined  as 

R0-i  (3.30) 

where  f  is  the  Coriolis  parameter  defined  as  20  sin  <f>,  0  is  the  rota¬ 
tional  speed  of  the  earth,  and  ^  is  the  latitude.  The  terms  inside  the 
parentheses  of  Equations  3.27-3,29  represent  nondimens ional  flows  (flows 
per  unit  depth  for  Equations  3.27  and  3.28).  The  terms  outside  the 
parentheses  (with  the  exception  of  AZ*)  convert  to  dimensional  flow. 

The  Rossby  scaling  arises  partially  from  the  fact  that  the  vector 
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potentials  contain  a  displacement  vector  computed  from  a  time  integral 
of  velocity. 

3.4  COMPUTATIONAL  SEQUENCE 

The  order  of  processor  computations  are  outlined  in  the  flowchart 
of  Figure  3.6  and  are  described  as  follows.  After  the  initializations 
within  CH3D,  Subroutine  WQMOUT  is  called  the  first  time  to  read  in  input 
parameters  and  the  map  file  required  by  the  processor.  Processor  ini¬ 
tializations  are  then  conducted,  and  time  invariant  geometric  informa¬ 
tion  (i.e.  layer  thickness,  horizontal  cell  lengths,  and  cell  surface 
areas)  are  computed  and  output.  Logic  returns  to  CH3D,  and  HM  time- 
stopping  commences. 

When  the  HM  time  iteration  counter,  IT,  equals  ITWQS  (the  iteration 
to  begin  processing  transport  information),  WQMCVOL  is  called.  This 
Call  statement  is  at  the  beginning  of  the  CH3D  time-stepi  ing  DO  loop 
(i.e.  the  beginning  of  the  HM  time  step).  WQMCVOL  is  an  entry  point 
within  WQMOUT  that  computes  and  outputs  cell  volumes  and  horizontal  flux 
facial  areas.  This  call  to  WQMCVOL  allows  the  HM  to  undergo  spin-up 
before  starting  transport  processing.  Volumes  and  facial  areas  at  that 
time  in  the  HM  simulation  are  required  as  initial  values  for  the  proc¬ 
essed  transport  information.  Logic  returns  to  CH3D,  and  hydrodynamic 
computations  are  performed  for  the  HM  time  step. 

If  IT  is  equal  to  or  greater  than  ITWQS,  WQTVD  is  called  at  the  end 
of  each  HM  computational  time  step  with  the  new  velocity  field.  WQTVD 
is  an  entry  point  within  WQMOUT  where  time-varying  transport  information 
is  processed. 
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CH3D 


Figure  3.6.  Interface  processor  flow  chart 
(Sheet  1  of  3) 
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Figure  3.6.  (Sheet  2  of  3) 


Figure  3.6.  (Sheet  3  of  3) 
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The  first  step  upon  entering  WQMOUT  at  WQTVD  is  to  update  time 
averages  of  velocities  and  the  cumulative  displacements  (i.e.  Equations 
3.20-3.22).  The  cumulative  displacements  are  time  integrations  of  the 
velocities  using  the  trapezoidal  rule.  Half  of  the  HM  time  step  is  used 
to  update  these  integrations  prior  to  updating  the  time  average  of  the 
cumulative  displacements  so  that  the  displacements  are  centered  within 
each  time  step  interval.  Subroutine  PROCZ  is  called  next  to  update  com¬ 
ponents  for  the  Stokes'  drifts  calculations. 

Subroutine  PROCZ  computes  spatial  averages  for  the  velocities  and 
the  cumulative  displacements  and  updates  the  time  averages  of  their  pro¬ 
ducts  (i.e.  the  vector  potentials  of  Equations  3.23-3.25).  Spatial  av¬ 
eraging  is  required  to  locate  the  velocities  and  cumulative  displace¬ 
ments  where  the  vector  potentials  are  defined  (see  Figure  3.5).  Logic 
is  then  returned  to  WQMOUT. 

Upon  returning  from  PROCZ,  the  second  half  of  the  HM  time  step  in¬ 
tegrations  for  the  cumulative  displacements  are  performed  to  advance 
these  integrations  to  the  end  of  the  HM  time  step  interval.  Following 
this  operation,  there  is  a  check  for  the  end  of  the  time -averaging  in¬ 
terval;  a  counter,  IKNT,  is  compared  with  the  NAVG,  the  number  of  HM 
time  steps  to  be  averaged.  If  IKNT  equals  NAVG,  Subroutine  PR0C2  is 
called;  otherwise  logic  proceeds  within  WQMOUT. 

PR0C2  is  an  entry  point  within  PROCZ  where  the  final  time -averaging 
updates  and  calculations  for  the  Stokes'  flow  components  are  performed. 
Spatial  averaging  is  performed  for:  velocities,  cumulative  displace¬ 
ments,  time-averaged  velocities,  time-averaged  cumulative  displacements, 
and  the  Jacobians  and  surface  layer  thicknesses  required  for  scaling  the 
A2  terms.  The  time -averaging  updates  for  the  vector  potentials  are 
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performed,  the  products  of  the  means  are  subtracted  (refer  to  Equation 
3.16),  and  the  results  are  scaled  by  the  Jacobians  as  in  Equation  3.26 
to  produce  the  A  terms.  Finally,  the  A  terms  are  set  to  zero  at  corners 
and  inflow  boundaries.  Logic  is  returned  to  WQMOUT  with  all  information 
required  to  compute  time -averaged  Stokes'  flows. 

Processor  computations  for  cell  flows  are  separated  for  horizontal 
and  vertical  directions.  WQMOUT  sweeps  through  all  horizontal  flow 
faces  checking  for  flow  direction  to  compute  the  horizontal  flow  with 
the  appropriate  velocity  and  flow  area.  The  time-averaged  horizontal 
flow  variable  (i.e.  the  Eulerian  residual)  is  updated,  and  the  horizon¬ 
tal  sweep  continues. 

After  computing  all  horizontal  flows ,  a  vertical  sweep  for  each 
column  is  made  to  compute  vertical  flows  and  vertical  diffusivities . 

The  time-averaged  vertical  flow  (i.e.  Eulerian  residual)  and  diffusivity 
are  also  updated  during  the  sweep. 

During  the  horizontal  and  vertical  sweeps  discussed  above,  IKNT  and 
NAVG  are  compared.  If  IKNT  and  NAVG  are  equal,  the  computations  for 
Stokes'  flows  are  completed  by  executing  Equations  3.27-3.29  in  finite 
difference  form  using  central  differences.  The  differences  A£  and  A f?  in 
Equations  3.27-3.29  are  unity  since  the  transformed  coordinates  are 
specified  as  integers  representing  coordinate  line  numbers  that  incre¬ 
ment  by  one. 

With  time-averaged  computations  complete,  all  time -averaged  vari¬ 
ables  (i.e.  horizontal  and  vertical  Eulerian  and  Stokes'  flows  and 
vertical  diffusivities)  are  written  to  an  output  file  for  later  use  in 
the  WQM.  The  cell  volumes  and  horizontal  flux  facial  areas  for  the  sur¬ 
face  boxes  ac  the  end  of  the  time -averaging  interval  are  also  computed 
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and  written  to  the  output  file  along  with  the  HM  time.  Finally,  all 
time -averaged  variables  are  reset  to  zero  for  preparation  of  the  next 
time -averaging  interval,  and  logic  is  returned  to  CH3D. 

It  is  noted  that  the  processor  is  independent  of  the  length  of  the 
averaging  interval,  thus,  intratidal  and  intertidal  averaging  can  be 
done  with  the  same  processor.  For  intratidal  averaging,  Stokes'  flows 
are  computed  and  output,  but  they  are  not  used  in  the  WQM  transport 
equation;  only  the  Eulerian  flows  are  used. 

3 . 5  CHAPTER  SUMMARY 

The  computational  methods  consist  of  the  HM,  WQM,  and  interface 
processor.  The  processor,  which  resides  within  the  HM,  translates  in¬ 
formation  from  the  HM  boundary- fitted  grid  to  the  integrated  compartment 
grid  of  the  WQM.  Time -invariant  geometric  information  and  time -varying 
transport  information  (i.e.  flows,  surface  cell  volumes,  and  vertical 
diffusivities)  are  computed,  converted  to  WQM  units,  and  output  while 
the  HM  is  executing.  This  output  file  is  subsequently  used  as  input  for 
driving  the  WQM  transport  equation. 

Computations  within  the  HM  are  nondimens ional  and  in  transformed 
coordinates.  Nondimens ional,  contravariant  velocities  are  converted  to 
dimensional,  physical  flows  for  WQM  use.  Conversion  from  contravariant 
to  physical  components  retains  flows  normal  to  transverse  grid  lines 
(i.e.  cell  faces)  as  required  by  the  WQM.  The  Stokes'  flows  are  com¬ 
puted  during  the  intertidal  averaging  interval,  thus  reducing  memory 
requirements.  Additionally,  all  averaging  for  displacements  and  veloci¬ 
ties  and  their  products  are  done  with  HM  variables  for  convenience,  then 
converted  to  Stokes'  flows  at  the  end  of  the  averaging  intervals.  The 
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processor  does  not  distinguish  from  intratidal  and  intertidal  averaging; 
rather,  the  WQM  uses  only  the  Eulerian  flows  for  intratidal  simulations 
and  adds  the  Stokes'  flows  to  the  Eulerian  flows  for  intertidal  simula¬ 
tions. 


CHAPTER  4 


EVALUATION  OF  RESIDUAL  COMPUTATIONS 

Application  of  Lagrangian  residual  computations  and  the  evaluation 
of  their  use  in  transport  simulations  are  presented  in  this  chapter. 
Results  are  presented  for  an  application  of  the  HM  and  interface  sub¬ 
routines  to  a  simple,  two-dimensional  (2D)  test  case  that  has  an  analyt¬ 
ical  result.  This  test  was  performed  to  verify  whether  the  methods  for 
computing  residual  currents  had  been  correctly  implemented.  The  HM  and 
WQM  were  next  applied  to  Chesapeake  Bay  to  evaluate  the  use  of  residual 
currents  for  salinity  transport.  The  results  of  this  application  and 
various  sensitivity  analyses  are  evaluated.  Finally,  a  discussion  of 
the  results  is  presented. 

4.1  ANALYTICAL  TEST  CASE 

The  2D  analytical  study  of  Ianniello  (1977)  was  chosen  to  test  the 
correctness  of  the  residual  computations.  Ianniello  (1977)  studied  the 
secondary  currents  generated  by  first-order  oscillatory  flows.  He  rec¬ 
ognized  that  the  steady-state  (residual)  secondary  currents  were  the 
difference  between  the  Lagrangian  and  Eulerian  residuals,  i.e.  the 
Stokes'  drift,  which  is  the  residual  velocity  caused  by  a  particle  pass¬ 
ing  through  spatially  vary-,,  *  velocity  fields  during  a  tidal  cycle. 

To  examine  tide  induced  residual  currents,  Ianniello  (1977)  analyt¬ 
ically  developed  solutions  of  the  residual  currents  (i.e.  Eulerian  resi¬ 
duals  and  Stokes'  drifts)  for  a  2D  (longitudinal  and  vertical), 
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rectangular  channel  of  constant  width  and  depth,  and  closed  at  one  end. 
The  solution  is  for  a  fluid  of  constant  density  and  no  wind  stress  on 
the  surface.  A  no  slip  condition  was  applied  to  the  bottom.  The  open 
end  was  forced  with  an  M2  tide,  and  a  no  flow  condition  was  imposed  at 
the  closed  end.  Therefore,  the  residual  currents  are  induced  solely  by 
the  nonlinear  interaction  of  the  tidal  currents. 

The  equations  solved  by  Ianniello  (1977)  were  the  nondimensional , 

2D  horizontal  momentum  and  continuity  equations  with  the  hydrostatic 
assumption  for  vertical  momentum.  The  solution  was  obtained  by  a  sec¬ 
ond-order  perturbation  analysis  for  weak  nonlinearity.  Ianniello  con¬ 
ducted  the  analysis  for  various  vertical  eddy  viscosity  models,  includ¬ 
ing  constant  eddy  viscosity.  Ianniello' s  solutions  provide  a  good  test 
case  for  checking  the  residual  currents  that  are  numerically  generated 
with  the  HM  and  the  interface  processor.  Similarly,  Najarian  et  al. 
(1982  and  1984)  used  Ianniello' s  results  to  test  their  residual  currents 
generated  with  a  2D,  laterally-averaged  hydrodynamic  model  (see  Section 
2.4). 

The  analytical  results  of  Ianniello  reported  by  Najarian  et  al. 
(1982  and  1984)  are  also  used  here.  These  results  were  obtained  with 
constant  vertical  eddy  viscosity,  a  linear  bottom  shear  stress  law,  >= 
1.0,  and  dQ  -  1.8,  where  and  d q  are  the  nondimensional  channel  length 
and  depth,  respectively.  For  a  tidal  period  of  12.42  hours  and  a  depth 
of  10  m,  1^  =  1.0  yields  a  dimensional  channel  length  of  about  70.0  km. 
The  value  for  the  eddy  viscosity  can  be  computed  from  the  rsl3tlcr«ship 
for  dQ  (Ianniello  1977);  dg  =  1.8  yields  a  vertical  eddy  viscosity,  Nz, 

O 

of  21.7  cnr/sec.  A  tidal  amplitude  of  30.0  cm  was  imposed.  A  HM  grid 
of  20  cells  long  by  10  cells  deep  was  used  to  represent  the  channel. 
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The  conditions  and  parameters  for  this  test  case  are  summarized  as 
follows : 

channel  length,  L  -  70.0  km 
channel  depth,  d  -  10  m 
AX  -  3.5  km 
AZ  -  1.0  ra 
At  -  276  sec 

tidal  amplitude  -  30.0  cm 
tidal  period  -  12.42  hours 
linear  bottom  shear  stress  law 

O 

vertical  eddy  viscosity,  Nz  -  21.7  cnr/sec 
no  horizontal  eddy  viscosity 

With  the  above  conditions,  the  HM  was  run  for  ten  tidal  cycles. 
Dynamic  steady-state  conditions  were  reached  after  about  two  or  three 
tidal  cycles.  Eulerian  residuals  and  Stokes'  drifts  were  computed  for 
the  tenth  tidal  cycle  and  were  saved  in  an  output  file.  Numerical  re¬ 
sults  were  compared  to  the  analytical  results  fur  the  vertical  profile 
of  residual  longitudinal  velocities  at  a  station  two  grid  cells,  or  7.0 
km,  from  the  open  boundary.  All  KM  runs  for  the  2D  test  case  were  con¬ 
ducted  on  the  WES  VAX  8800. 

The  first  runs  of  the  HM  were  made  with  the  vertical  length  scale, 
Zr  (also  referred  to  as  ZREF) ,  set  to  1000  cm.  This  parameter  is  used 
by  the  HM  for  nondimensionalization  (see  Appendix  A).  A  value  of  1000 
cm  is  appropriate  to  provide  a  nonaimensionai  depth  on  the  order  of  1.0 
as  recommended  by  model  documentation  (Sheng  1986b) .  With  a  value  of  Zr 
-  1000  and  Nz  -  21.7  cm^/sec,  the  numerically  generated  residual  cur¬ 
rents  did  not  closely  match  the  analytical  result  as  shown  in  Figure 
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4.1.  Experimentation  with  the  HM  revealed  that  adjustment  of  the  param¬ 
eter  Zr  modified  the  computed  residual  currents.  With  a  value  of  Zr  ~ 
500,  close  agreement  with  the  analytical  result  was  obtained  as  shown  in 
Figure  4.2.  Further  decreases  in  Zr  required  increasing  the  value  of  Nz 
above  the  analytical  value  to  obtain  satisfactory  agreement.  For  ex¬ 
ample,  with  a  value  of  Zr  -  250,  Nz  had  to  be  increased  to  a  value  of 
around  30.0  cm^/sec  to  provide  satisfactory  results  as  shown  in  Figure 
4.3.  To  obtain  good  agreement  with  Zr  -  1000,  the  value  of  Nz  had  to  be 
decreased  considerably.  Figure  4.4  presents  results  for  Zr  -  1000  and 
Nz  -  13.0. 

The  dependence  of  the  numerical  results  on  the  value  of  Zr  was 
first  thought  to  be  due  to  an  error  in  the  HM  code  related  to  the  non- 
dimensional  relationships.  However,  after  several  independent  searches 
for  coding  errors,  no  coding  errors  were  found.  It  is  speculated  that 
the  dependence  of  the  HM  on  Zr  is  a  result  of  numerical  diffusion.  The 
test  results  shown  in  Figures  4. 1-4.4  indicate  that  as  Zr  is  increased, 
Nz  must  be  decreased  to  align  numerical  results  with  analytical  results. 
Numerical  diffusion  (i.e.  artificial  dissipation  or  dampening)  is  pro¬ 
portional  to  the  product  of  velocity  and  the  spatial  step  size  for  up¬ 
wind  differencing  schemes  (Roach  1972)  such  as  those  used  for  the  con¬ 
vective  acceleration  terms  of  the  HM.  Numerical  vertical  diffusion  of 
momentum  is  proportional  to  the  dimensional  vertical  velocity  and  layer 
thickness,  which  suggests  a  Zr^  dependence  (refer  to  Appendix  A).  Nu¬ 
merical  dissipation  has  the  same  effect  as  increasing  the  real  eddy  vis¬ 
cosity  or  diffusiv'.ty.  Therefore,  Nz  had  to  be  reduced  from  the  analyt¬ 
ical  value  when  using  the  recommended  value  of  Zr  -  10 
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VERTICAL  EDDY  VISCOSITY  =  21.70  sq  cm/sec  . . „  ANALYTICAL  RESULT 
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LAGRANGIAN  RESIDUAL,  cm/s 


Figure  4.1.  Comparison  of  numerical  model  with  analytical  result, 

2D  tidal  channel  test,  N  «  21.7,  Z  **  1000 
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Figure  4.2.  Comparison  of  numerical  model  with  analytical  result, 
2D  tidal  channel  test,  N  =  2~.7,  Z  *»  500 
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2D  ANALYTICAL  TEST  CASE 

VERTICAL  EDDY  VISCOSITY  =  30.00  sq  cm/sec 
ZREF  =  250.00 
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Figure  4.3.  Comparison  of  numerical  model  with  analytical  result, 

2D  tidal  channel  test,  N  =  30.0,  Z  =  250 
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2D  ANALYTICAL  TEST  CASE 

VERTICAL  EDDY  VISCOSITY  =  13.00  sq  cm/sec 
ZREF  =  1000.00 


ANALYTICAL  RESULT 
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Figure  4.4.  Comparison  of  numerical  model  with  analytical  result, 
2D  tidal  channel  test,  N  =  13.0,  Z  *»  1009 
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It  should  be  noted  in  Figures  4. 1-4. 4  that  the  residual  velocities 
of  the  numerical  model  were  not  plotted  for  the  surface  layer  because 
the  values  were  off  the  scale.  The  large  values,  in  opposite  direc¬ 
tions,  for  the  Eulerian  residuals  and  Stokes'  drifts  are  caused  by  the 
fully  conservative  condition  imposed  in  the  surface  layer.  Recall  from 
Chapter  3  that  a  completely  conservative  formulation  is  used  to  compute 
the  residual  currents.  To  guarantee  conservation  of  flow  and  mass  in 
the  surface  layer,  there  can  not  be  any  flux  through  the  water  surface, 
which  means  the  velocities  and  vector  potentials  must  be  zero  on  the 
surface.  Mass  conserving  Eulerian  flows  in  the  surface  layer  are  guar¬ 
anteed  by  averaging  flows  (i.e.  velocity  times  an  area  that  accounts  for 
water  surface  displacement),  rather  than  velocities.  The  analytical 
model  was  not  concerned  with  a  mass  transport  application,  and  a  conser¬ 
vative  formulation  for  the  surface  layer  was  not  imposed.  Therefore, 
the  two  results  differ  in  the  surface  layer  for  the  Eulerian  and  Stokes 
components.  However,  when  the  two  numerical  components  are  added,  the 
resulting  Lagrangian  current  shows  close  agreement  with  the  analytical 
result. 

A  non- conservative  formulation  in  the  processor  was  tested  by  not 
forcing  conservative  Eulerian  residuals  in  the  surface  layer  (i.e.  simp¬ 
ly  average  the  velocities)  and  computing  vector  potentials  on  the  water 
surface.  The  latter  was  accomplished  by  assuming  a  zero  vertical  gradi¬ 
ent  of  horizontal  velocity  at  the  surface  and  using  the  surface  dis¬ 
placement  over  the  HM  time  step  as  the  vertical  velocity  on  the  surface. 
This  test  yielded  surface  layer  Eulerian  and  Stokes  components  close  to 
the  analytical  result. 
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The  Lagrangian  velocity  profiles  of  Figures  4. 1-4.4  sum  to  zero 
from  top  to  bottom.  This  result  is  consistent  with  the  fact  that  there 
is  no  net  change  in  volume  of  the  channel,  and  there  is  no  inflow  at  the 
closed  end.  The  tides  induce  a  residual  circulation  in  the  channel,  but 
inflow  and  outflow  should  balance  resulting  in  no  net  flow  into  or  out 
the  channel.  Both  the  analytical  and  numerical  residual  flows  do  bal¬ 
ance.  It  is  interesting  to  note  that  if  just  the  Eulerian  residual  cur¬ 
rents  are  considered,  there  would  appear  to  be  a  net  outflow,  which  is 
impossible.  Positive  velocities  are  out  of  the  channel  towards  the 
ocean. 

Although  the  residual  currents  for  this  test  case  are  quite  small 
compared  with  the  intratidal  currents,  which  peak  on  the  order  of  30.0 
cm/sec,  it  is  important  to  recognize  that  the  Stokes'  drift  is  on  the 
same  order  as  the  Eulerian  residual  currents.  These  results  well  il¬ 
lustrate  the  need  to  consider  Lagrangian  residuals  rather  than  Eulerian 
residuals  when  concerned  with  net  transport.  With  the  relatively  close 
agreement  between  the  aodel/processor  results  and  the  analytical  re¬ 
sults,  the  computational  procedures  of  Chapter  3  were  concluded  to  be 
properly  implemented. 


4.2  CHESAPEAKE  BAY  APPLICATION 

The  computational  system,  consisting  of  the  HM,  the  interface  pro¬ 
cessor  within  the  HM,  and  the  WQM,  was  applied  to  Chesapeake  Lay  to 
evaluate  its  usefulness  for  real  transport  applications.  Chesapeake  Bay 


is  a  partially  mixed  estuary  of  the  drowned  river  valley  type  (i.e. 
coastal  plain).  The  Bay  is  approximately  300  km  long  and  50  km  wide  at 
the  widest  point.  The  average  depth  is  about  eight  meters.  Freshwater 
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inflows  are  delivered  by  eight  major  rivers  of  which  the  James,  Potomac, 
and  Susquehanna  Rivers  provide  about  82  percent  of  the  flow  (HydroQual 
1987) .  The  tidal  range  is  the  greatest  at  the  mouth  of  the  Bay  with  an 
average  of  about  0.8m.  Tidal  currents  average  0.5  knots  with  peaks  as 
high  as  2.0  knots  (HydroQual  1987).  Winds  can  have  a  significant  effect 
on  tidal  heights  and  currents. 

The  Chesapeake  Bay  application  is  divided  into  three  parts:  Septem¬ 
ber  1983  Simulation,  1985  Simulation,  and  Sensitivity.  September  1983 
was  one  of  the  periods  used  for  HM  calibration;  the  year  1985,  a  rela¬ 
tively  dry  year,  was  used  for  HM  verification. 

4.2.1  September  1983  Simulation 

The  September  1983  simulation  provided  a  verification  of  the  HM  and 
WQM  interfacing  procedures  for  Chesapeake  Bay.  As  part  of  the  Chesa¬ 
peake  Bay  model  study  (Dortch  et  al.  1988  and  Dortch  1988),  the  HM  was 
calibrated  for  September  1983  because  of  the  availability  of  tide,  ve¬ 
locity,  and  salinity  data.  The  results  of  the  HM  calibration  have  been 
reported  by  Kim  et  al.  (1989)  and  show  excellent  agreement  with  observa¬ 
tions.  September  1983  was  a  particularly  interesting  test  case  because 
of  a  strong  wind  mixing  event  that  reduced  salinity  stratification 
followed  by  re-stratification  near  the  end  of  the  month.  September  1983 
was  considered  a  severe  test  case  of  the  HM  and  the  interfacing  proce¬ 
dures.  This  relatively  short  simulation  period  also  provided  a  practi¬ 
cal  means  of  conducting  initial  model  interfacing  tests  without  exces¬ 
sive  disk  storage  requirements  and  computation  time. 

4. 2. 1.1  Intratidal  Tests.  The  interface  processor  was  first  test¬ 
ed  through  intratidal  transport  comparisons  of  HM  and  WQM  salinity. 
Transport  comparisons  were  facilitated  by  the  fact  that  both  the  HM  and 


74 


WQM  can  transport  a  conservative  dissolved  substance  (e.g.  salinity). 

The  tests  consisted  of  comparing  salinity  computed  with  the  WQM  against 
salinity  computed  with  the  HM.  The  WQM  should  not  be  any  more  accurate 
that  the  HM  when  comparing  with  observed  salinity  since  the  WQM  uses  HM 
information  for  transport.  In  fact,  the  HM  does  not  have  to  be  cali¬ 
brated  in  order  to  use  it  as  a  standard  for  testing  the  ability  of  the 
WQM  to  accurately  use  HM  information.  Observed  salinity  data  were 
omitted  from  the  salinity  plots  presented  in  this  section  to  facilitate 
comparison  of  WQM  against  HM  salinity. 

For  intratidal  transport,  the  WQM  should  yield  transport  results 
that  are  practically  identical  to  that  of  the  HM,  assuming  the  models 
are  properly  coupled,  the  solutions  schemes  are  similar,  and  the  in¬ 
tratidal  averaging  interval  of  HM  information  is  not  so  large  that  re¬ 
sidual  currents  become  important  (i.e.  Stokes'  drift).  Intratidal 
transport  tests  included  use  of  HM  information  without  any  time  averag¬ 
ing  and  time  averaging  of  HM  information  over  one  and  two  hour  periods. 

The  model  grid  used  for  Chesapeake  Bay  is  shown  in  Figure  3.1. 

There  are  729  surface  cells.  The  number  of  layers  vary  from  a  minimum 
of  two  in  the  shallows  to  a  maximum  of  15  in  the  deep  trench,  resulting 
in  a  total  of  3948  computational  cells .  The  WQM  grid  is  identical  to 
the  HM  grid  except  that  the  HM  has  five  additional  columns  along  the 
ocean  boundary  since  the  WQM  grid  at  the  ocean  boundary  must  start  one 
row  inside  the  HM  boundary  (see  Section  3.3). 

The  calibrated  HM  and  its  data  for  September  1983  were  furnished  by 
the  WES  Hydraulics  Laboratory.  To  test  the  model  interfacing,  the  HM 
had  to  be  re-run  with  the  processor  included.  The  HM  was  started  on 
September  1,  1983  with  observed  temperature  and  salinity  distributions. 
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The  model  was  run  with  the  temperature  and  salinity  fields  held  constant 
for  five  days  to  allow  the  velocity  field  to  develop  from  initially  zero 
to  conditions  reflecting  the  water  density  variations.  After  the  five 
day  spin-up  period,  the  salinity  and  temperature  fields  were  released, 
and  the  model  was  run  for  another  21  days  with  the  changing  salinity  and 
temperature  affecting  the  hydrodynamics  and  vice  versa  (i.e.  baroclinic 
coupling) .  Monthly  averaged  boundary  conditions  were  used  for  river 
inflow  rates  and  temperatures  and  ocean  boundary  salinity.  Ocean 
boundary  salinity  did  vary  across  the  width  and  depth  of  the  boundary. 
Time -varying  observed  tides  were  used  at  the  ocean  boundary.  The  HM  and 
processor  were  executed  on  a  Cray  2  located  at  the  U.S.  Array  Tank  Auto¬ 
motive  Command  in  Warren,  Michigan,  via  a  satellite  link.  The  26  day  HM 
simulation  required  about  half  an  hour  of  CPU  time  on  the  Cray  2  with  a 
ten  minute  time  step. 

The  processing  of  time-varying  transport  information  was  started  at 
the  beginning  of  HM  simulation  day  6  and  was  continued  until  the  end  of 
the  HM  run.  Salinity  data  computed  by  the  HM  were  stored  on  disk  at  30 
minute  intervals  as  the  HM  was  executing,  and  transport  information  for 
the  WQM  was  stored  at  the  end  of  each  averaging  interval.  Subsequently, 
the  WQM  was  executed  on  the  Cray  2  for  the  21  day  simulation  period 
using  the  previously  computed  and  stored  transport  information. 

Salinity  data  computed  by  the  WQM  were  stored  on  disk  at  a  frequen¬ 
cy  equal  to  the  WQM  time  step.  The  salinity  data  computed  by  the  HM  and 
WQM  were  plotted  together  for  comparison  versus  time  for  three  layers 
(surface,  mid-depth,  bottom)  at  multiple  stations.  Stations  used  for 
salinity  comparisons  are  shown  in  Figure  4.5.  The  station  noted  PR  in 
the  September  1983  simulation  plots  is  Station  LP  in  Figure  4.5. 
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The  WQM  salinity  values  were  practically  identical  to  HM  salinity 
when  HK  information  was  processed  without  any  averaging  (i.e.  flow  up¬ 
dates  for  transport  and  model  time  steps  were  the  same,  ten  minutes). 
This  test  indicated  that  the  two  models  were  successfully  interfaced, 
and  the  differences  in  solution  schemes  did  not  cause  any  differences  in 
salinity  transport  for  the  21  day  simulation.  At  the  time  of  this  work, 
the  HM  did  not  have  the  QUICKEST  solution  scheme  for  salinity  advection. 
Therefore,  the  upwind  differencing  option  was  used  in  the  WQM  for  hori¬ 
zontal  advection  for  all  of  the  September  1983  tests. 

Salinity  computed  by  the  HM  (CH3D)  and  the  WQM  is  compared  at  se¬ 
lected  stations  in  Figure  4.6  for  an  intratidal  averaging  period  of  one 
hcur  (i.e.  averaging  over  six  HM  time  steps).  The  WQM  time  step  for 
this  simulation  was  also  fixed  at  one  hour.  Only  Eulerian  residuals 
were  used  for  the  intratidal  WQM  simulations.  Examination  of  Figure  4.6 
reveals  that  the  WQM  salinity  closely  tracked  the  HM  salinity.  A  slight 
phase  shift  on  the  same  order  as  the  averaging  interval  is  observable  in 
Figure  4.6.  The  WQM  results  were  not  adjusted  to  account  for  the  phase 
lag  created  by  the  averaging  interval.  The  greatest  difference  in  WQM 
and  HM  salinity  occurred  at  the  Bay  entrance  (Station  0C1  and  0C2)  just 
after  150  hours.  The  WQM  tended  to  overshoot  the  HM  result  briefly. 

This  was  due  to  a  short-lived  stability  violation  that  quickly  damped 
out  without  numerically  diverging.  The  one  hour  time  step  of  the  WQM 
temporarily  exceeded  the  horizontal  flow  stability  criterion. 

The  HM  was  also  run  with  the  processor  set  for  two  hour  averaging, 
and  the  WQM  was  run  with  this  input  and  a  one  hour  time  step.  The  re¬ 
sults  were  very  similar  to  those  shown  in  Figure  4.6,  except  for  an  in¬ 
crease  in  the  phase  lag  of  the  averaging  interval.  However,  infci'atJ.  lal 
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Figure  4.6,  Salinity  computed  with  intratidal  HM  and  intratidal 
WQM  for  September  1983  (Sheet  1  of  8) 


CONCENTRATION,  PPT  CONCENTRATION,  PPT  CONCENTRATION,  PPT 


79 


STATION  OC2 


SURFACE 


MID-DEPTH 


30 
25 
20 
15 
10 
5 
0 

100  150  200  250  300  350  400  450  500  550  500  65G 

TIME,  HRS 


BOTTOM 

v  u  v-  ,vl'V\/VVVV  'i^JvUv  w\a/'Vv:  Vv-  vN1  ''-V-v 


i—i— L-j-i  .1  i _ i _ I- 1  i _ 1  I  i-i-.l  ■ !  »  I  i  -i-.  I _ i _ 1—1— l-jl~L  i—i 


Figure  4.6.  (Sheet  2  of  8) 
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Figure  4.6 
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Figure  4.6.  (Sheet  6  of  8) 
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Figure  4.6.  (Sheet  7  of  8) 
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Figure  4.6.  (Sheet  8  of  8) 
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averaging  over  three  hours  produced  WQM  salinity  that  did  not  track  HM 
salinity  nearly  as  well  as  the  runs  with  shorter  averaging  periods. 
Therefore,  intratidal  averaging  (i.e.  where  only  Eulerian  residuals  are 
used)  should  be  limited  to  about  two  hours  or  less  for  semi-diurnal 
tides.  Because  of  the  nearly  identical  agreement  of  the  HM  and  WQM 
intratidal  salinity  transport  results ,  it  was  concluded  that  the  WQM  was 
properly  interfaced  with  the  HM. 

4. 2. 1.2  Intertidal  Test.  The  September  1983  simulation  was  re¬ 
peated  with  an  intertidal  averaging  period  of  12.5  hours  (i.e.  the  aver¬ 
aging  period  contained  75  HM  time  steps).  A  constant  WQM  time  step  of 
3750  seconds,  which  divides  evenly  into  12.5  hours,  was  used  for  the  WQM 
run.  It  was  possible  to  use  a  larger  time  step  for  the  intertidal  run 
than  for  the  intratidal  run  because  the  residual  flows  are  smaller  for 
averaging  periods  extending  over  the  full  tidal  period.  It  would  have 
been  possible  to  use  an  even  larger  WQM  time  step  had  it  not  been  for 
the  period  of  high  winds  for  sustained  periods  during  the  middle  of  the 
month.  The  Lagrangian  residuals  (Eulerian  residuals  plus  Stokes 'drifts) 
were  used  for  the  intertidal  WQM  run. 

The  comparison  of  WQM  and  HM  salinity  for  the  intertidal  transport 
test  is  presented  in  Figure  4.7.  From  these  plots,  it  is  evident  that 
the  tides  have  been  averaged  out  of  the  WQM  solution.  The  WQM  salinity 
have  been  obtained  through  more  gradually  varying  residual  currents. 
Therefore,  the  WQM  should  not  yield  the  high  frequency  tidally  varying 
salinity  that  the  KM  does,  but  should  generally  follow  the  mean  of  the 
HM  salinity. 

At  the  Bay  entrance,  the  WQM  salinity  tends  to  be  slightly  shifted 
near  the  troughs  of  the  HM  salinity  at  the  southern  end  (Station  0C1) , 
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Figure  4.7.  Salinity  computed  with  intratidal  HM  and  intertidal 
WQM  for  September  1983  (Sheet  1  of  8) 
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Figure  4.7.  (Sheet  2  of  8) 
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Figure  4.7.  (Sheet  3  of  8) 
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Figure  4.7.  (Sheet  4  of  8) 
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Figure  4.7.  (Sheet  5  of  8) 
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Figure  4.7.  (Sheet  6  of  8) 
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Figure  4.7.  (Sheet  7  of  8) 
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especially  for  the  mid-depth  and  bottom  layers.  However,  at  the  north¬ 
ern  end  of  the  entrance  (Station  0C3),  the  WQM  salinity  tends  to  be 
closer  to  the  peaks  of  the  HM  salinity,  especially  for  the  surface  and 
mid-depth.  These  trends  suggest  that  there  is  a  net  flow  out  (yielding 
lower  salinity)  at  the  southern  end,  and  there  is  a  net  flow  into  the 
Bay  at  the  Northern  end  (yielding  higher  salinity) .  This  observation  is 
in  general  agreement  with  observed  similar  salinity  skewness  in  the  Bay 
near  the  mouth.  Such  net  flows  and  resulting  skewness  can  be  attributed 
to  the  Coriolis  force  which  creates  a  counterclockwise  circulation  that 
pushes  down-estuary  flow  towards  the  western  shore,  thus  exiting  at  the 
southern  end  of  the  mouth.  Likewise,  up-estuary  flow  is  deflected  tow¬ 
ard  the  eastern  shore  and  would  tend  to  enter  at  the  northern  end  of  the 
mouth . 

The  WQM  salinity  closely  follows  the  trends  of  the  HM  salinity  in 
Figure  4.7.  The  results  do  not  compare  as  closely  as  do  the  intratidal 
results  because  of  the  loss  of  some  detail  of  the  hydrodynamic  informa¬ 
tion  when  averaged  over  an  intertidal  period.  Recall  from  Chapter  2 
that  the  first-order  approximation  for  Lagrangian  residuals  does  not 
account  for  higher-order  effects  (i.e.  tidal  phase  dispersion).  The  WQM 
does  follow  the  HM  even  during  the  intense  wind  mixing  event  (around 
hour  500)  followed  by  re-stratification  (see  the  mid-depth  plot  at  Bay 
Bridge,  Station  BB  of  Figure  4.7). 

The  results  of  the  21  day  intertidal  simulation  were  encouraging, 
but  the  simulation  period  was  not  long  enough  to  adequately  test  the 
success  of  the  methodology.  Chesapeake  Bay  is  so  large  relative  to  the 
inflows  and  tidal  flows,  that  it  can  take  months  for  the  effects  of  re¬ 
sidual  currents  to  have  a  substantial  impact  on  salinity  distributions. 
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A  year' long  simulation  provides  the  real  test  of  transport  computations 
using  residual  currents. 

4.2.2  1985  Simulation 

During  the  Chesapeake  Bay  model  study,  the  HM  was  verified  by  the 
WES  Hydraulics  Laboratory  against  observed  salinity  data  over  the  entire 
year  1985.  Salinity  data,  which  had  been  collected  at  approximately  two 
week  intervals  during  data  gathering  cruises  by  various  local  univer¬ 
sities,  were  extracted  from  the  Chesapeake  Bay  Program  data  base  (US 
Environmental  Protection  Agency  1989).  The  HM  and  input  data  sets  for 
1985  were  furnished  by  the  WES  Hydraulics  Laboratory  following  model 
verification. 

Daily  updates  were  used  for  flow  and  temperature  at  river  bound¬ 
aries  throughout  1985.  River  flows  were  obtained  from  recorded  gages, 
and  equilibrium  temperatures  (Edinger  et  al.  1974)  were  computed  from 
meteorological  data  and  used  for  river  inflow  temperatures.  Continuous 
(i.e.  one  hour  interval)  tidal  elevation  records  at  the  Bay  entrance 
were  used  to  drive  the  tidal  boundary.  Observed  salinity  data  for  the 
Bay  entrance  at  approximately  two  week  intervals  were  used  for  the  ocean 
boundary;  values  between  observations  were  linearly  interpolated.  Sal¬ 
inity  data  observed  throughout  the  Bay  in  early  January  1985  were  inter¬ 
polated  over  the  HM  grid  and  used  for  i.nitial  conditions  to  spin-up 
(i.e.  start-up)  the  model. 

WES  brought  on-line  an  in-house  Cray  Y/MP  supercomputer  during  this 
phase  of  the  study.  A  HM  simulation  of  the  full  year  1985  required  ap¬ 
proximately  six  CPU  hours  on  the  Cray  Y/MP.  The  HM  required  a  five 
minute  time  step  for  1985  to  maintain  stability  throughout  the  year.  HM 
information  was  averaged  over  12.5  hour  periods  (i.e.  approximately  a 
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tidal  period  which  contained  150  HM  time  steps)  and  was  subsequently 
used  for  the  WQM  intertidal  simulation.  At  the  time  that  the  1985  simu¬ 
lations  were  conducted,  the  QUICKEST  scheme  for  advection  of  salinity 
had  been  implemented  in  the  HM  in  all  three  directions.  The  QUICKEST 
differencing  was  also  used  for  horizontal  advection  in  the  WQM;  however, 
fully  implicit,  central  differencing  (Euler  implicit  method)  was  used 
for  vertical  advection. 

The  autostepping  (i.e.  automatic  time  stepping)  feature  was  used 
for  the  WQM  1985  simulations.  Although  tidally- averaged  flow  updates 
were  used,  time  steps  less  than  a  tidal  cycle  were  required  to  maintain 
stability  for  the  explicit  horizontal  advection  scheme.  Therefore,  the 
WQM  used  a  constant  field  of  hydrodynamics  for  each  WQM  time  step  until 
it  was  time  for  the  next  hydrodynamic  field  update.  The  average  WQM 
time  step  for  the  1985  simulation  was  3,182  sec  with  a  maximum  allowable 
time  step  of  8,693  sec.  These  are  an  order  of  magnitude  greater  than 
the  300  sec  time  step  used  by  the  HM.  The  smallest  WQM  time  step  was  3 
sec,  which  was  the  time  step  required  to  take  the  WQM  time  exactly  to 
the  next  update  interval  for  rime- averaged  HM  input.  The  total  CPU  time 
required  by  the  WQM  for  the  1985  salinity  simulation  was  488  seconds  on 
the  Cray  Y/MP. 

Observed,  HM,  and  WQM  salinity  were  compared  for  the  1985  simula¬ 
tion.  Comparisons  are  presented  with  the  same  format  used  in  Section 
4.2.1,  i.e.  salinity  versus  time  for  surface,  mid- depth,  and  bottom  at 


utuj.  tipic 


a  uotiwita  . 


Seven  more  stations  were  added  for  a  more  thorough 


comparison.  Salinity  comparisons  for  the  intratidal  HM,  intertidal  WQM 
and  observations  are  presented  in  Figure  4.8  for  all  stations  shown  in 
Figure  4.5.  Observed  salinity  data  were  not  available  for  all  depths 
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Figure  4.8.  Salinity  computed  with  intratidal  HM  and  intertidal 
WQM  (base  conditions)  for  1985  (Sheet  1  of  15) 


CONCENTRATION,  PPT  CONCENTRATION,  PPT  CONCENTRATION,  PPT 


99 


LEGEND 


STATION  OC2  (CB  7.4) 
SURFACE 


.  WQM 

_  CH3D 

•  OBSERVED 


35 
30 
25 
20 
15 
10 
5 
0 

0  50  100  150  200  250  300  350  400 

TIME,  DAYS 


MID-DEPTH 


BOTTOM 


Figure  4.8.  (Sheet  2  of  15) 
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4.8.  (Sheet  3  of  15) 


CONCENTRATION,  PPT  CONCENTRATION.  PPT  CONCENTRATION.  PPT 


101 


LEGEND 

....  WQM 

_  CH3D 

•  OBSERVED 


35 
30 
25 
20 
15 
10 
5 
0 

0  50  100  150  200  250  300  350  400 

TIME,  DAYS 

MID-DEPTH 

35 

30 
25 
20 
15 
10 
5 
0 

0  50  100  150  200  250  300  350  400 

TIME,  DAYS 

BOTTOM 

35 
30 
25 
20 
15 
10 
5 
0 

0  50  100  150  200  250  300  350  400 

TIME,  DAYS 


STATION  WT  (CB  6.3) 
SURFACE 


Figure  4.8 


(Sheet  4  of  15) 


CONCENTRATION,  PPT  CONCENTRATION.  PPT  CONCENTRATION,  PPT 


102 


STATION  MB  (CB  5.1) 
SURFACE 


UE&EhLP 

.  WQM 

_  CH3D 

•  OBSERVED 


MID-DEPTH 


BOTTOM 


Figure  4.8.  (Sheet  5  of  15) 
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Figure  4.8.  (Sheet  6  of  15) 
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Figure  4.8.  (Sheet  7  of  15) 
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Figure  4.8.  (Sheet  8  of  15) 
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and  stations  presented.  Salinity  data  computed  by  the  HM  and  WQM  were 
averaged  over  12.5  hour  periods  for  plotting  to  facilitate  comparisons. 
Otherwise,  the  intratidal  HM  salinity  fluctuations  are  large  enough  to 
obscure  WQM  results  and  observations  when  plotted  on  the  scale  of  Figure 
4.8  for  an  entire  year.  WQM  salinity  data  were  also  averaged  over  12.5 
hour  periods  to  be  consistent  with  HM  results.  Stations  that  contain 
only  two  plots  were  two  layers  deep.  The  numbers  in  parentheses  next  to 
the  station  identification  are  additional  location  information  that 
either  corresponds  to  the  WQM  surface  box  number  (e.g.  B  71)  or  a  Chesa¬ 
peake  Bay  observation  station  (e.g.  CB  6.3). 

Examination  of  the  plots  of  Figure  4.8  reveals  the  ability  of  the 
intertidal  WQM,  using  residual  currents,  to  track  the  HM  salinity.  The 
WQM  results  contain  more  variability  than  those  obtained  from  the  HM, 
even  after  tidally  averaging  the  computed  salinity.  This  is  attributed 
to  the  use  of  tidally  averaged  HM  information  in  the  WQM;  flows  may 
change  substantially  from  one  hydrodynamic  update  to  the  next.  Flows  in 
the  HM  changed  more  gradually  (i.e.  five  minute  updates),  thus,  result¬ 
ing  in  more  gradual  changes  in  salinity. 

Although  the  WQM  does  not  reproduce  every  detail  of  the  HM  salini¬ 
ty,  it  does  follow  the  trends  over  the  entire  year  quite  well.  Some  of 
the  salinity  values  change  significantly  (e.g,  10  part  per  thousand, 
ppt)  over  a  period  of  a  few  days  as  shown  at  the  Bay  Bridge  (Station  BB 
of  Figure  4.8).  The  ability  of  the  WQM  to  follow  the  general  transport 
character  of  the  HM  at  the  upper  Bay  station  (Station  UB  of  Figure  4.8) 
is  also  quite  impressive. 

Several  statistics  were  computed  to  assist  in  quantitatively  evalu¬ 
ating  how  well  the  WQM  tracked  the  HM  for  various  simulation  results. 
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For  all  42  locations  (fifteen  stations  with  two  to  three  layers  each) 
where  results  are  presented,  the  differences  in  WQM  and  HM  salinity 
(i.e.  residual  of  WQM  minus  HM,  or  model  error)  were  obtained.  The  mean 
error  (ME) ,  the  mean  absolute  error  (MAE) ,  the  root  mean  square  of  the 
errors  (RMSE) ,  and  the  standard  deviation  (SD)  of  the  MAE  were  calcu¬ 
lated  for  all  stations  combined.  ME  indicates  the  sign  of  error,  MAE 
gives  the  average  error  without  regard  to  sign,  and  RMSE  is  the  standard 
deviation  of  error  about  zero  mean  error.  The  SD  of  the  MAE  was  com¬ 
puted  so  that  hypotheses  concerning  two  means  could  be  tested  for  sig¬ 
nificance  (Miller  and  Freund  1977).  All  significance  testing  reported 
herein  used  a  -  0.01  level,  and  the  sample  sizes  were  corrected  for 
autocorrelation  (Reckhow  et  al.  1986).  Statistics  for  the  simulation 
results  of  Figure  4.8  were  -0.84,  1.21,  and  1.76  ppt  for  ME,  MAE,  and 
RMSE,  respectively.  All  reported  residual  statistics  are  summarized  in 
Table  4.1. 

The  tributaries  contain  the  greatest  fresh  and  salt  water  gradi¬ 
ents,  thus,  they  are  the  most  difficult  areas  of  the  Bay  to  model.  The 
greatest  deviations  of  computed  HM  salinity  from  observed  salinity  occur 
in  the  tributaries.  With  the  relatively  coarse  grid  scale  employed 
here,  it  is  not  surprising  that  the  models  do  not  accurately  match  ob¬ 
served  data.  The  WQM  salinity  also  differ  the  most  from  those  computed 
by  the  HM  in  the  tributaries.  Differences  in  transport  properties  (i.e. 
intratidal  versus  intertidal),  model  time  step  size,  and  the  vertical 


advection  solution  schemes  could  lead  to 


differences  in  the  two  models. 


This  issue  is  discussed  in  more  detail  in  the  Discussion  section. 

The  models  match  observed  salinity  well  in  the  upper  and  lower  Po¬ 
tomac  River  (Stations  UP  and  LP).  However,  neither  model  matches 
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Table  4.1.  Summary  of  residual  statistics 


Test  Condition 

ME 

MAE 

RMSE 

£D 

Intertidal  WQM  base  condi¬ 
tions  (Fig.  4.8) 

-0.84 

1.21 

1.76 

1.26 

Without  Stokes'  flows 
(Fig.  4.9) 

-3.02 

3.13 

4.50 

1.97 

Without  vertical  diffusion 
(Fig.  4.10) 

0.62 

1.75 

2.82 

2.04 

Hydrodynamics  averaged  over 
25  hrs  (Fig.  4.12) 

-1.15 

1.50 

2.10 

1.45 

With  horizontal  diffusion  of 
100  m2/s  (Fig.  4.13) 

-0.99 

1.28 

1.88 

1.15 

WQM  time  step  -  300  sec 

-0.82 

1.23 

1.81 

1.30 

Intratidal  WQM  base  condi¬ 
tions  (Fig.  4.17) 

-0.48 

0.74 

1.28 

0.86 

Intertidal  vs  intratidal 

-0.31 

0.98 

1.51 

1.30 

WQM  (Fig.  4.18) 


Note:  ME  -  Mean  residual  or  error 
MAE  -  Mean  absolute  error 
RMSE  -  Root  mean  square  error 
SD  -  Standard  deviation  of  MAE 
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observed  salinity  very  well  in  mid-Potomac  River  (Station  MP  of  Figure 
4.8).  This  disparity  is  attributed  to  inadequate  vertical  grid  resolu¬ 
tion.  The  mid- Potomac  River  reach  is  typically  seven  to  ten  meters 
deep,  whereas,  only  two  layers  (about  four  meters)  were  used  in  most  of 
the  reach.  Therefore,  the  observed  bottom  salinity  data  for  this  reach 
were  measured  in  water  much  deeper  than  the  model's  bottom.  The  HM  and 
WQM  should  have  included  more  layers  for  the  Potomac  River. 

4.2.3  Sensitivity 

4. 2. 3.1  Stokes'  Flows.  The  first  sensitivity  test  consisted  of 
running  the  intertidal  WQM  without  the  Stokes'  flows  (i.e.  using  only 
Eulerian  residual  flows) .  The  same  conditions  for  the  1985  simulation 
were  repeated  except  that  the  Stokes'  flows  were  not  added  to  the  Euler¬ 
ian  residuals  for  the  advective  terms  during  WQM  execution. 

The  results  of  the  simulation  without  Stokes'  flows  are  presented 
in  Figure  4.9.  After  about  a  month  of  simulation,  it  is  apparent  that 
the  WQM  does  not  satisfactorily  reproduce  the  HM  results  or  observed 
results  at  most  stations.  The  only  stations  where  the  WQM  agrees  well 
with  observations  and  HM  results  are  Stations  UP  and  UJ,  where  fresh 
water  flows  persist.  The  results  indicate  that  a  considerable  amount  of 
salt  water  is  missing  without  the  Stokes'  flows.  Therefore,  the  effect 
of  the  Stokes'  flows  is  to  transport  more  salinity  into  the  estuary. 

The  WQM  versus  HM  error  statistics  for  this  simulation  were  -3.02,  3.13, 
and  4.50  for  ME,  MAE,  and  RMSE,  respectively.  These  statistical  results 
are  somewhat  misleading.  Main  bay  errors,  which  were  smaller  those  of 
the  tributaries,  tend  to  draw  down  the  relatively  large  error  statistics 
of  the  tributaries  and  upper  bay.  The  results  of  Figures  4.8  and  4.9 
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clearly  demonstrate  the  need  to  properly  compute  residual  currents  for 
intertidal  transport  modeling. 

4. 2. 3. 2  Vertical  Diffusion.  A  simulation  was  made  without  any 
vertical  diffusion  in  the  WQM  to  test  the  importance  of  this  mechanism 
in  the  intertidal  transport  model.  Tidally  averaged  vertical  diffusiv- 
ities  ranged  from  a  minimum  of  0.005  cm^/sec  to  a  maximum  of  about  1000 
cmVsec*  Again  the  same  1985  conditions  were  used,  and  the  WQM  was 
driven  with  Lagrangian  residual  currents  (i.e.  Eulerian  residuals  and 
Stokes'  flows).  Rather  than  using  the  time-averaged  vertical  diffusiv- 
ities  output  by  the  processor,  the  vertical  diffusivity  was  set  to  zero 
in  the  WQM.  The  results  of  this  run,  which  are  shown  in  Figure  4.10, 
indicate  that  the  WQM  does  not  track  observed  data  or  HM  results  nearly 
as  well  as  those  with  the  tidally-averaged  vertical  diffusivities  (Fig¬ 
ure  4.8).  In  general,  the  WQM  over-predicted  salinity  in  the  bottom 
layers.  The  HM-WQM  error  statistics  for  this  run  were  0.62,  1.75,  and 
2.82  for  the  ME,  MAE,  and  RMSE,  respectively.  Although  the  disparity  is 
not  as  great  as  that  obtained  by  neglecting  the  Stokes'  flows  (Figure 
4.9),  the  results  do  indicate  the  need  to  include  vertical  diffusivities 
computed  by  the  HM. 

4. 2. 3. 3  Upwind  Differencing,  The  importance  of  using  a  high-order 
of  accuracy  for  the  horizontal  advection  was  tested  by  making  HM  and  WQM 
runs  with  the  first-order  upwind  differencing  scheme.  The  HM  and  WQM 
results  of  Figure  4.8  were  obtained  with  QUICKEST,  whereas  the  HM  and 
WQM  results  of  Figure  4.11  were  obtained  with  upwind  differencing  for 
horizontal  advection  of  salinity.  For  vertical  advection,  the  QUICKEST 
scheme  was  still  used  in  the  HM,  and  the  Euler  implicit  method  was  used 
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Figure  4.11.  Salinity  computed  with  intratidal  HM  and  intertidal  WQM 
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The  greatest  differences  of  the  upwind  results  (i.e.  Figure  4.11) 
from  the  QUICKEST  results  (i.e.  Figure  4.8)  occurred  in  the  upper  bay 
and  upper  tributaries  (Stations  UB,  UJ,  UP,  and  UY) ,  where  the  HM  over¬ 
predicts  salinity  with  the  low-order,  upwind  scheme.  The  upwind  scheme 
is  numerically  diffusive  and  diffuses  salinity  upstream.  Upwind  dif¬ 
ferencing  tends  to  inhibit  the  ability  of  the  WQM  to  track  the  HM  sali¬ 
nity  even  in  the  main  bay  over  a  long  simulation  period.  Differences  in 
HM  vs  WQM  results  is  attributed  to  differing  amounts  of  numerical  dam¬ 
pening  (i.e.  diffusion)  induced  by  the  two  models.  Error  statistics 
were  not  used  for  this  test  since  the  test  was  run  to  examine  the  effect 
of  low- order  differencing  for  horizontal  advection  rather  than  examine 
HM-WQM  residuals  due  to  low-order  differencing. 

From  the  upwind  differencing  results,  it  is  apparent  that  the  high¬ 
er-order  accurate  scheme  is  more  important  for  the  HM  than  for  the  WQM. 
This  is  a  reasonable  assessment  since  numerical  dampening  of  the  upwind 
scheme  is  proportional  to  the  product  of  velocity  and  the  spatial  step 
size.  The  spatial  resolution  is  the  same  for  the  HM  and  WQM,  but  the 
intertidal  velocities  of  the  WQM  are  considerably  less  than  the  in- 
tratidal  velocities  of  the  HM.  Therefore,  numerical  dampening  of  the 
upwind  scheme  is  greater  in  the  HM  than  in  the  WQM.  With  the  relatively 
coarse  grid  used  here,  it  is  not  possible  to  resolve  regions  with  large 
salinity  gradients  (i.e.  upper  bay  and  tributaries)  in  the  HM  without  a 
high-order  advection  scheme. 

Although  salinity  is  modeled  reasonably  well  with  an  upwind  scheme 
in  the  intertidal  WQM,  other  water  quality  constituents  with  stronger 
concentration  gradients  could  exhibit  more  than  a  desirable  amount  of 
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numerical  diffusion.  It  is  prudent  to  use  the  relatively  inexpensive, 
but  effective,  higher-order  QUICKEST  scheme. 

4. 2. 3. 4  Averaging  Interval.  All  intertidal  transport  tests  dis¬ 
cussed  above  used  a  12.5  hour  averaging  interval,  which  is  approximately 
the  average  tidal  period  of  Chesapeake  Bay.  It  is  possible  that  longer 
averaging  periods,  such  as  several  tidal  periods,  could  still  capture 
the  transport  characteristics  of  the  mean  flows  and  tide-induced  flows. 
An  averaging  period  of  25.0  hours  (i.e.  approximately  two  tidal  cycles) 
was  tested  to  investigate  the  effect  of  longer  averaging  time.  The  HM 
was  run  with  NAVG  set  to  300  iterations,  i.e.  the  processor  averaged 
over  300  HM  time  steps.  The  WQM  was  subsequently  run  with  these  inter¬ 
tidal  Hydrodynamics.  The  results  of  this  simulation  are  shown  in  Figure 
4.12.  Comparison  with  Figure  4.8  indicates  little  loss  of  the  proper 
salinity  transport  trends  with  the  longer  averaging  interval.  However, 
error  statistics  of  ME  -  -1.15,  MAE  -  1.50,  and  EMSE  -  2.10  are  greater 
than  those  associated  with  the  single  tidal  period  averaging  interval. 
The  two  means  (i.e.  MAEs  for  single  and  two  tidal  cycle  averaging)  are 
significantly  different  (a  *»  0.01). 

If  the  time -averaging  interval  is  extended  too  far,  the  intertidal 
transport  could  fail  to  resolve  changes  arising  from  the  non- tidal  forc¬ 
ing,  such  as  major  shifts  in  winds  which  occur  at  about  three  day  inter¬ 
vals,  sea  state  changes,  and  time-varying  fresh  water  inflows.  For 
time-varying  applications,  ;s  suggested  that  the  averaging  interval 
should  not  exceed  two  or  three  tidal  cycles.  The  benefits  gained  (i.e. 
smaller  processed  files)  with  longer  averagir*"  periods  are  not  justified 
due  to  the  loss  in  transport  resolution.  In  all  cases,  intertidal 
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Figure  4.12.  Salinity  computed  with  intratidal  HM  and  intertidal  WQM 
(with  hydrodynamics  averaged  over  25.0  hrs)  for  1985 
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transport  results  should  be  compared  with  intratidal  transport  to  ensure 
that  the  proper  transport  characteristics  have  been  captured. 

4. 2. 3. 5  Horizontal  Diffusion.  All  HM  and  WQM  simulations  dis¬ 
cussed  above  did  not  have  any  physical  horizontal  diffusion.  Physical 
horizontal  diffusion  was  considered  to  be  orders  of  magnitude  smaller 
than  advective  transport  of  salinity  and  was  not  necessary  in  HM  salini¬ 
ty  calibration.  Additionally,  diffusion  associated  with  spatial  averag¬ 
ing  (i.e.  dispersion)  is  less  important  in  3D  models.  A  run  was  made 
where  physical  horizontal  diffusion  was  set  to  100.0  m  /sec  in  the  WQM 
to  test  the  sensitivity.  It  is  important  to  recognize  that  this  value 
of  horizontal  diffusion  is  estimated  to  be  one  to  two  orders  of  mag¬ 
nitude  less  than  that  introduced  by  the  upwind  differencing  scheme  for 
advection.  Horizontal  diffusion  in  HM  was  still  zero.  Selected  sta¬ 
tions  for  the  run  with  horizontal  diffusion  are  shown  in  Figure  4.13. 
Results  with  horizontal  diffusion  appeared  to  match  the  HM  in  the  main 
bay  about  as  well  as  the  results  without  diffusion  shown  in  Figure  4.8; 
this  observation  is  exhibited  by  the  salinity  plot  of  Station  MB  (Figure 
4.13)  which  is  very  similar  to  the  results  without  diffusion  (i.e.  Fig¬ 
ure  4.8).  In  the  upper  bay  (Station  UB)  and  in  the  upper  tributaries 
(e.g.  Station  UY) ,  diffusion  tends  to  cause  the  WQM  to  more  closely 
match  the  HM  as  shown  in  Figure  4.13.  However,  in  the  middle  to  lower 
reaches  of  the  tributaries,  the  WQM  under-predicted  HM  salinity  more 
with  diffusion  (see  Station  MJ  of  Figure  4,13).  The  WQM  diffusion  was 
applied  evenly  throughout  the  grid;  thus,  diffusion  increases  salinity 
in  the  raid- tributaries  and  decreases  salinity  in  the  lower  tributaries. 

The  WQM  vs  HM  residual  statistics  for  the  simulation  with  horizon¬ 
tal  diffusion  were  -0.99,  1.28,  and  1.88  for  the  ME,  MAE,  and  RMSE, 
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Figure  4.13.  Salinity  computed  with  intratidal  HM  and  intertidal  WQM 
(with  horizontal  diffusion  set  to  100  m^/sec)  for  1985 
(Sheet  1  of  4) 
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respectively.  The  WQM  results  with  rather  large  horizontal  diffusion 
are  graphically  similar  to  the  base  case  (i.e.  no  diffusion).  However, 

the  difference  in  the  MAEs  for  this  simulation  and  the  simulation  with¬ 

out  diffusion  are  significant,  and  it  can  be  concluded  that  the  WQM 
results  with  diffusion  did  not  match  the  HM  as  well  as  the  results  with¬ 
out  diffusion. 

4. 2. 3. 6  Time  Step  Size.  The  sensitivity  of  the  WQM  to  the  time 
step  size  was  tested.  The  WQM  simulation  of  1985  was  conducted  with  the 

WQM  time  step  set  300  sec,  the  same  value  used  for  the  HM.  The  plotted 

sal*~ity  appeared  nearly  identical  to  the  WQM  salinity  of  Figure  4.8, 
which  were  computed  with  the  much  larger  time  steps  that  occurred  with 
autostepping.  The  MAE  for  this  run  was  1.23,  which  is  not  significantly 
different  (a  =  0.01)  from  the  MAE  of  the  results  in  Figure  4.8. 

The  result  that  the  WQM  is  insensitive  to  time  step  size  is  rea¬ 
sonable  since  the  QUICKEST  scheme  used  for  horizontal  advection  is 
third-order  accurate  in  space  and  time.  Although  the  Euler  implicit 
method  used  for  WQM  vertical  advection  is  only  first-order  accurate  in 
time,  the  numerical  dampening  that  is  introduced  is  quite  small.  This 
is  due  to  the  fact  that  the  dampening  is  proportional  to  the  product  of 
the  velocity  squared  and  the  time  step  size  (Anderson  1984).  The  verti¬ 
cal  velocity  is  very  small  (i.e.  on  the  order  of  mm/sec,  or  smaller), 
thus,  the  dampening  is  small  for  the  time  steps  utilized. 

4.3  DISCUSSION 

4.3.1  Characteristics  of  Residual  Currents 

4. 3. 1.1  2D  Test  with  Salinity  and  River  Flow.  The  1985  simula¬ 
tions  with  and  without  the  Stokes'  flows  operable  indicated  that  the 
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Stokes'  flows  tend  to  transport  more  salt  water  into  the  Bay.  The  2D 
analytical  test  case  of  Section  4.1  showed  that  the  Stokes'  drift  ad- 
vects  water  up  the  estuary  away  from  the  mouth.  However,  the  idealized 
results  of  Section  4.1  were  obtained  without  salinity  variations  or 
fresh  water  inflow  as  occur  in  Chesapeake  Bay  and  most  other  real  es¬ 
tuaries.  Several  2D  simulations  for  the  idealized  estuary  of  Section 
4.1  were  re-run  with  conditions  of  varying  salinity  and  fresh  water 
(i.e.  river)  inflow  to  determine  if  the  Stokes'  drift  still  had  the  ef¬ 
fect  of  forcing  water  up-estuaryt  thus,  pushing  salinity  into  the  es¬ 
tuary.  It  is  well  known  that  salinity  stratification  (i.e.  density 
stratification)  coupled  with  fresh  water  inflow  induces  gravitational 
circulation  with  surface  water  moving  down- estuary  toward  the  mouth  and 
bottom  water  moving  up-estuary  (Dyer  1973  and  Officer  1976).  However, 
it  is  not  obvious  what  the  direction  of  the  Stokes'  drift  is  under  these 
conditions . 

The  2D  tests  were  run  with  salinity,  with  and  without  river  inflow. 
These  runs  were  carried  out  20  tidal  cycles  instead  of  10  to  achieve 
dynamic  steady- state.  The  ocean  salinity  boundary  condition  was  held 
constant  to  the  initial  values  set  at  the  mouth.  Values  for  vertical 
eddy  viscosity  and  ZREF  were  set  to  those  obtained  during  calibration 
with  the  analytical  result  of  Section  4.1.  However,  it  is  recognized 
that  density  stratification  suppresses  vertical  mixing.  The  results 
shown  in  Figure  4.14  were  obtained  without  river  inflow  and  with  initial 
salinity  values  of  20.0  ppt  on  the  surface  increasing  linearly  to  29.0 
ppt  on  the  bottom;  there  was  not  any  horizontal  variation  in  salinity 
initially.  Upon  reaching  dynamic  steady- state,  there  was  a  slight  ver¬ 
tical  stratification  of  about  one  half  ppt.  Results  are  compared  with 
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2D  ANALYTICAL  TEST  CASE  WITH  SALINITY 
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Figure  4.14.  Numerical  model  result  for  2D  tidal  channel  test  with 

salinity,  without  river  inflow,  N  -  13.0 
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the  analytical  result  of  no  salinity  and  no  river  inflow.  From  Figure 
4.14,  it  is  evident  that  the  Stokes'  drift  is  still  in  the  up-estuary 
direction.  The  same  test  was  also  run  with  a  smaller  value  for  the  ver- 

O 

tical  eddy  viscosity  of  1.0  cm  /sec.  The  smaller  viscosity  had  the  ef¬ 
fect  of  increasing  the  steady-state  vertical  stratification  (one  to  two 
ppt  variation  vertically),  which  resulted  in  tide-induced  gravitational 
circulation  as  evident  in  Figure  4.15.  There  is  still  an  up-estuary 
Stokes'  flow  on  the  bottom.  There  is  also  an  up-estuary  tide-induced 
Eulerian  residual  flow  on  the  bottom.  It  is  interesting  that  tides 
alone  can  induce  gravitational  residual  circulation.  The  vertical  sum 
of  the  Lagrangian  residuals  is  zero  for  both  Figures  4.14  and  4.15  as  it 
should  be  without  river  inflow. 

O 

The  2D  test  was  next  run  with  a  river  inflow  of  about  143  nr /sec 
(5000  ft3/sec) .  Salinity  was  initialized  with  a  linear  variation  of  1 
ppt  to  20  ppt  horizontally  with  each  column  vertically  increasing  by  1 
ppt  per  layer;  thus,  the  ocean  varied  from  20  ppt  on  the  surface  to  30 
ppt  on  the  bottom.  Horizontal  salinity  variations  were  used  to  reduce 
the  time  to  reach  dynamic  steady- state,  at  which  point  about  one  half 
ppt  vertical  stratification  was  remaining.  The  typical  gravitational 
circulation  of  a  partially  mixed  estuary  is  evident  in  Figure  4.16  for 
both  the  Eulerian  and  Lagrangian  residuals.  The  vertical  sum  of  the 
Lagrangian  residuals  is  not  zero,  rather  it  corresponds  to  the  river 
inflow.  The  up-estuary  Stokes'  drift  is  still  present  with  the  river 
inflow.  It  is  now  clear  that  tide-induced  Stokes'  drift  tends  to  advect 
water  up-estuary  regardless  of  whether  or  not  there  are  baroclinic  forc¬ 
ing  (density- induced  flows)  and  river  inflow.  Therefore,  neglecting 
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2D  ANALYTICAL  TEST  CASE  WITH  SALINITY 
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Figure  4.15.  Numerical  model  result  for  2D  tidal  channel  test  with 
salinity,  without  river  inflow,  Nz  =  1.0 
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Figure  4.16.  Numerical  model  result  for  2D  tidal  channel  test  with 
salinity,  with  river  inflow,  N  -  13.0 
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Stokes'  flows  in  intertidal  transport  would  tend  to  under-estimate  salt 
water  intrusion,  which  is  illustrated  by  Figure  4.9. 

4. 3. 1.2  Residual  Flows  and  Velocities.  Computed  residual  flows 
and  velocities  were  examined  at  several  of  the  Chesapeake  Bay  stations. 
The  1985  intertidal  hydrodynamic  data  output  by  the  processor  (i.e.  hor¬ 
izontal  and  vertical  Eulerian  residuals  and  Stokes'  flows)  were  averaged 
over  the  year-long  simulation  period  to  examine  the  average  magnitudes 
of  various  components.  Stations  MB,  UB,  and  MJ  were  selected  for  this 
examination.  Annual  averages  of  the  tidally  averaged  currents  (i.e. 
volumetric  flow  rate  and  velocity),  in  both  the  horizontal  (i.e.  along 
the  primary  flow  axis)  and  vertical  directions  for  the  Eulerian, 

Stokes',  and  Lagrangian  residuals,  are  presented  in  Tables  4. 2-4.4  for 
Stations  MB,  UB,  and  MJ,  respectively.  Annual  averages  of  the  tidally 
averaged  hydrodynamic  information  provide  summaries  of  the  relative  mag¬ 
nitudes  of  the  currents  without  the  noise  of  the  daily  fluctuations  as¬ 
sociated  with  shifting  winds  and  changing  flows  and  tides. 

Examination  of  Tables  4. 2-4.4  reveals  several  interesting  features 
concerning  the  nature  of  the  residual  currents.  The  direction  of  flow 
at  all  three  stations  is  similar  to  the  baroclinic  flow  patterns  ob¬ 
tained  in  the  2D  tests  (i.e.  Figures  4.14-4.16),  where  horizontal  Euler¬ 
ian  and  Lagrangian  residuals  are  down-estuary  on  the  surface  and  up-es- 
tuary  on  the  bottom.  Additionally,  the  direction  of  the  horizontal 
Stokes'  drift  tends  to  transport  salinity  up-estuary.  The  velocities 
are  of  the  same  order  of  magnitude  as  those  obtained  with  the  2D  test. 

Although  the  Stokes'  drift  is  in  general  smaller  than  the  Eulerian 
residuals,  it  is  apparent  that  the  tide- induced  Stokes'  flows  are  not 
negligible  for  long-term  transport  computations.  The  horizontal  Stokes' 
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Table  4.2.  1985  annual  averages  of  intertidal  flows  and  velocities 

at  Station  MB 

Horizontal 


O 

Flowfm  /sec)  Velocity  (’em/sec) 


Laver  Eulerian 

Stokes 

Lagrangian 

Eulerian 

Stokes 

Lagrangian 

1  (surface) 

-545.1 

-55.5 

-600.6 

-8.38 

-0.81 

-9.19 

2 

-353.2 

-12.5 

-365.7 

-8.86 

-0.31 

-9.17 

3 

-242.1 

23.3 

-218.8 

-5.07 

0.49 

-4.58 

4 

-145.3 

23.0 

-122.3 

-3.04 

0.48 

-2.56 

5 

-47.0 

22.8 

-24.2 

-0.98 

0.48 

-0.50 

6 

61.4 

24.8 

86.2 

1.29 

0.52 

1.81 

7 

154.4 

30.8 

185.2 

3.23 

0.64 

3.87 

8 

225.7 

32.7 

258.4 

4.73 

0.68 

5.41 

9 

275.9 

30.4 

306.3 

5.78 

0.64 

6.42 

10 

298.3 

26.0 

324.3 

6.25 

0.54 

6.79 

11 

287.5 

19.0 

306.5 

6.02 

0.40 

6.42 

12 

251.2 

13.6 

264.8 

5.26 

0.28 

5.54 

13 

228.7 

10.8 

239.5 

4.79 

0.23 

5.02 

14 

200.3 

9.4 

209.7 

4.19 

0.20 

4.39 

15  (bottom) 

160.0 

6.7 

166.7 

3.35 

0.14 

3.49 

Note:  Flows  and  velocities  are  in  the  £  (North-South)  direction. 
Positive  is  up-estuary,  negative  is  down-estuary. 

Vertical 

Flowfm /sec)  Velocity (mm/sec) 


Laver  Eulerian 

Stokes 

Lagrangian 

Eulerian  Stokes 

Lagrangian 

1  (surface) 

-16.5 

-56.5 

-73.0 

-0.0015  -0.0051 

-0.0066 

2 

-10.5 

-83.4 

-93.9 

-0.0009  -0.0075 

-0.0084 

3 

8.9 

-101.9 

-93.0 

0.0008  -0.0091 

-0.0083 

4 

38.4 

-113.4 

-75.0 

0.0034  -0.0102 

-0.0068 

5 

70.6 

-118.6 

-48.0 

0.0063  -0.0106 

-0.0043 

6 

94.4 

-109.6 

-15.2 

0.0085  -0.0098 

-0.0013 

7 

114.0 

-94.9 

19.1 

0.0102  -0.0085 

0.0017 

8 

130.0 

-76.7 

53.3 

0.0116  -0.0069 

0.0047 

9 

139.4 

-55.1 

84.3 

0.0125  -0.0049 

0.0076 

10 

139.5 

-30.4 

109.1 

0.0125  -0.0027 

0.0098 

11 

129.3 

-7.7 

121.6 

0.0116  -0.0007 

0.0109 

12 

85.2 

-3.4 

81.8 

0.0076  -0.0003 

0,0073 

13 

51.2 

-2.4 

48.8 

0.0046  -0.0002 

0.0044 

14  (bottom) 

22.7 

-1.6 

21.1 

0.0020  -0.0001 

0.0019 

Note:  Positive  is  up,  toward  surface. 
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Table  4.3.  1985  annual  averages  of  intertidal  flows  and  velocities 

at  Station  UB 


Horizontal 

FlowCrn^/sec")  Velocity  (’cm/sec’> 


Laver  Eulerian 

Stokes 

Lagraneian 

Eulerian 

Stokes 

Lagranzian 

1  (surface) 

-552.4 

-110.1 

-662.5 

-12.15 

-2.29 

-14.44 

2 

-360.6 

-19.7 

-380.3 

-11.08 

-0.60 

-11.68 

3 

-268.0 

-22.5 

-290.5 

-8.24 

-0.68 

-8.92 

4 

-150.4 

-22.5 

-172.9 

-4.62 

-0.69 

-5.31 

5 

-28.5 

3.1 

-25.4 

-0.87 

0.10 

-0.77 

6 

69.2 

17.2 

86.4 

2.13 

0.53 

2.66 

7  (bottom) 

86.3 

16.3 

102.6 

2.65 

0.50 

3.15 

Note:  Flows  and  velocities  are  in  the  £  (North- South)  direction. 
Positive  is  up-estuary,  negative  is  down-estuary. 


Vertical 


Flowfm^/sec)  Velocitv(nun/sec') 


Laver  Eulerian 

Stokes 

Laeraneian 

Eulerian  Stokes 

Lagraneian 

1  (surface) 

176.8 

-130.0 

46.8 

0.0210  -0.0155 

0.0055 

2 

250.6 

-205.8 

44.8 

0.0298  -0.0245 

0.0053 

3 

262.3 

-248.3 

14.0 

0.0312  -0.0295 

0.0017 

4 

220.2 

-223.5 

-3.3 

0.0262  -0.0266 

-0.0004 

5 

129.8 

-128.4 

1.4 

0.0154  -0.0153 

0.0001 

6  (bottom) 

57.3 

-47.3 

10.0 

0.0068  -0.0056 

0.0012 

Note:  Positive  is  up,  toward  surface. 
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Table  4.4.  1985  annual  averages  of  intertidal  flows  and  velocities 

at  Station  MJ 


Horizontal 


O 

Flowfnr /sec)  Velocitvf cm/sec) 


Laver  Eulerian 

Stokes 

Laeraneian 

Eulerian 

Stokes 

Laeraneian 

1  (surface) 

-275.5 

-43.0 

-318.5 

-5.75 

-1.06 

-6.81 

2 

-177.4 

121.3 

-56.1 

-5.51 

3.77 

-1.74 

3 

-7.3 

89.2 

81.9 

-0.23 

2.77 

2.54 

4  (bottom) 

74.7 

27.0 

101.7 

2.32 

0.84 

3.16 

Note:  Flows  and  velocities  are  in  the  rj  (East-V?est)  direction. 
Positive  is  up-estuary,  negative  is  down-estuary. 


Vertical 


FlowCm^/sec) 

Velocity (mm/sec) 

Laver  Eulerian 

Stokes 

Laeraneian 

Eulerian  Stokes 

Laeraneian 

1  (surface) 

-48.5 

-51.6 

-100.1 

-0.0043  -0.0046 

-0.0089 

2 

77.0 

-19.5 

57.5 

0.0069  -0.0017 

0.0052 

3  (bottom) 

31.0 

-2.5 

28.5 

0.0028  -0.0002 

0.0026 

Note:  Positive  is  up,  toward  surface. 
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flows  are  an  order  of  magnitude  smaller  than  the  Eulerian  flows  at  Sta¬ 
tion  MB,  but  they  are  of  about  a  half  to  the  same  order  of  magnitude  at 
Stations  UB  and  MJ. 

The  Lagrangian  horizontal  flows  at  Station  MJ  sum  to  -191.0  m^/sec, 
which  is  of  the  same  order  of  magnitude  as  the  annual  average  flow  in 

O 

the  James  River  for  1985,  about  325  m  /sec.  Recall  from  the  2D  tests 
that  the  net  Lagrangian  flow  (i.e.  the  sum  over  the  vertical)  in  a  tidal 
channel  is  the  freshwater  inflow.  The  model  has  three  cells  across  the 
channel  at  Station  M J ,  so  some  of  the  James  River  flow  passes  through 
the  adjacent  cells. 

Although  the  vertical  velocities  are  two  to  three  orders  of  mag¬ 
nitude  smaller  than  the  horizontal  velocities,  the  vertical  Eulerian  and 
Stokes'  flows  are,  in  general,  about  the  same  order  of  magnitude  as  the 
horizontal  flows.  The  vertical  Eulerian  residuals  are  generally  upward, 
while  the  vertical  Stokes'  drift  is  downward;  the  two  components  are 
about  the  same  order  of  magnitude,  thus  having  a  cancelling  effect. 

4.3.2  Salinity  Transport  Comparisons 

The  intertidal  WQM  predicted  salinity  lower  than  the  HM  in  the  tri¬ 
butaries  (refer  to  Figure  4.8).  Some  differences  in  WQM  and  HM  salini¬ 
ty  may  arise  from  the  fact  that  the  approach  for  obtaining  Lagrangian 
residual  currents  from  the  HM  is  only  a  first-order  approximation.  Re¬ 
call  from  Chapter  2  that  for  regions  with  considerable  nonlinearities 
(i.e.  where  k  >  0.5),  second-order  tidal  phase  dispersion  may  become 
significant  (Gomez-Reyes  1989).  Values  for  k  (based  on  mean  tidal  am¬ 
plitude  and  model  depth)  were  evaluated  for  various  portions  of  the  tri¬ 
butaries.  In  the  upper  reaches  of  the  Potomac,  James,  and  York  Rivers, 
where  tidal  amplitude  is  relatively  large  (compared  with  lower  reaches 
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of  the  tributaries)  and  depth  is  shallow,  k  is  in  the  range  of  0.12  to 
0.16.  The  values  of  k  in  the  upper  tributaries  are  considerably  higher 
than  for  the  whole  Bay  average  (i.e.  k  =  0.05),  but  they  are  consider¬ 
ably  lower  than  0.50,  the  value  that  Gomez-Reyes  (1989)  suggested  as  the 
acceptable  upper  bound  for  using  the  first-order  estimates  for  Lagran- 
gian  residuals.  However,  Gome:  Reyes  defined  k  as  the  ratio  of  tidal 
excursion  to  the  length  scale  of  the  velocity  gradient  for  flow  around  a 
headland,  which  may  have  little  meaning  in  this  study. 

There  was  concern  that  differences  in  WQM  and  HM  salinity  could 
indicate  a  failing  of  the  first-order  intertidal  averaging  method  to 
fully  resolve  transport.  To  address  this  question,  intratidal  hydrody¬ 
namics  (i.e.  averaging  over  15  HM  time  steps,  or  4500  sec)  were  pro¬ 
cessed  for  1985,  and  the  WQM  was  run  in  the  intratidal  mode  without  the 
need  for  Stokes'  flows.  The  results  of  this  simulation  (shown  in  Figure 
4.17)  when  compared  with  Figure  4.8  indicate  that  intratidal  averaging 
produces  tighter  agreement  with  the  HM  for  most  stations,  especially  in 
the  main  bay  and  lower  tributaries.  The  HM-WQM  error  statistics  for  the 
intratidal  run  were  -0.48,  0.74,  and  1.28  for  the  ME,  MAE,  and  RMSE,  re¬ 
spectively.  The  MAE  for  this  run  is  significantly  smaller  than  the  MAE 
for  the  results  in  Figure  4.8.  Therefore,  tidal  averaging  does  result 
in  some  loss  in  transport  information,  probably  second-order  tidal  phase 
effects  that  are  not  accounted  for  in  the  first-order  estimates. 

There  are  some  differences  in  the  HM  and  WQM  salinity  in  the  mid¬ 
dle  and  upper  tributaries  and  in  the  upper  bay  for  the  intratidal  simu¬ 
lation  of  1985  (see  Stations  UB,  MJ,  MP,  and  UY  of  Figure  4.17).  In¬ 
tratidal  salinity  computed  with  the  WQM  should  compar ;  closely  with 
those  computed  by  the  HM  at  all  stations  unless  the  differences  in  the 
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Figure  4.17.  Salinity  computed  with  intratidal  HM  and  intratidal 
WQM  (base  conditions)  for  1985  (Sheet  1  of  15) 
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Figure  4.17.  (Sheet  2  of  15) 
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Figure  4.17.  (Sheet  4  of  15) 
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Figure  4.17.  (Sheet  5  of  15) 
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Figure  4.17.  (Sheet  6  of  15) 
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Figure  4.17.  (Sheet  7  of  15) 
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Figure  4.17.  (Sheet  10  of  15) 
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Figure  4.17.  (Sheet  11  of  15) 
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Figure  4.17.  (Sheet  14  of  15) 
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solution  schemes  impact  the  results.  Differences  in  salinity  of  Figure 
4.17  could  indicate  that  the  intertidal  WQM  and  HM  salinity  differences 
for  corresponding  areas  of  Figure  4.8  may  not  totally  be  due  to  loss  of 
transport  information  during  the  averaging  procedure.  The  intertidal 
transport  results  at  Stations  UB,  MJ,  MP,  and  UY  compare  about  as  well 
with  the  HM  as  do  the  intratidal  results . 

The  major  difference  in  the  solution  schemes  is  that  the  HM  used 
QUICKEST  for  vertical  advection,  whereas  the  WQM  used  the  Euler  implicit 
method.  This  potential  source  of  model  disparity  was  investigated  by 
plotting  intertidal  WQM  versus  intratidal  WQM  salinity.  Any  differ¬ 
ences  in  solution  schemes  are  immediately  removed,  and  the  differences 
in  results  are  due  solely  to  the  length  of  the  averaging  interval.  The 
intertidal  WQM  salinity  of  Figure  4.8  and  the  intratidal  WQM  salinity  of 
Figure  4.17  are  plotted  together  in  Figure  4.18.  Comparison  of  Figures 
4.17  and  4.18  for  Stations  MJ,  LP,  MP,  and  UY  illustrates  the  improve¬ 
ments  realized  by  removing  differences  in  model  solution  schemes.  Error 
(i.e.  intertidal  minus  intratidal  salinity)  statistics  of  Figure  4.18 
were  -0.31,  0.98,  and  1.51  for  the  ME,  MAE,  and  RMSE,  respectively.  The 
only  source  of  model  disparity  in  Figure  4.18  is  the  effect  of  tidal 
averaging. 

There  are  two  primary  sources  of  model  disparity  in  Figure  4.8, 
tidal  averaging  of  the  hydrodynamics  and  differences  in  HM  and  WQM  solu¬ 
tion  schemes.  The  primary  source  of  model  disparity  in  Figure  4.17  is 
the  difference  in  the  HM  and  WQM  solution  schemes.  The  MAEs  were  1.21 
and  0.74  for  results  of  Figures  4.8  and  4.17,  respectively,  and  all 
three  MAEs  (i.e.  for  Figures  4.8,  4.17,  and  4.18)  are  significantly  dif¬ 
ferent  (a  »  0.01).  Although  the  disparity  (i.e.  MAE)  caused  by  tidal 
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Figure  4.18.  Salinity  computed  with  intratidal  WQM  and  intertidal 
WQM  for  1985,  base  conditions  (Sheet  1  of  15) 
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Figure  4.18.  (Sheet  3  of  15) 
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Figure  4.18.  (Sheet  4  of  15) 
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Figure  4.18.  (Sheet  5  of  15) 
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Figure  4.18.  (Sheet  6  of  16) 
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Figure  4.18.  (Sheet  7  of  15) 
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Figure  4.18.  (Sheet  8  of  15) 
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Figure  4.18.  (Sheet  10  of  15) 
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Figure  4.18.  (Sheet  11  of  15) 
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Figure  4.18.  (Sheet  14  of  15) 
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averaging  is  greater  than  that  caused  by  differences  in  HM  and  WQM  solu¬ 
tion  schemes,  the  latter  definitely  contributes  to  the  MAE  of  Figure  4.8 
results. 

Based  on  the  above  observations,  it  can  be  concluded  that  some  of 
the  differences  in  the  WQM  and  HM  salinity,  especially  in  the  middle  and 
upper  tributaries,  are  a  result  of  differences  in  the  numerical  solution 
schemes  of  the  two  models.  It  can  also  be  concluded  that  first-order, 
Lagrangian  residual  transport  differs  some  from  intratidal  transport. 
Intertidal  transport  introduced  an  average  absolute  error  of  less  than 
1.0  ppt  for  the  15  stations  compared.  These  differences  are  most  likely 
a  result  of  tidal  phase  dispersion,  which  is  not  accounted  for  in  the 
first-order  approach. 

The  root  cause  of  the  solution  scheme  disparity  associated  with 
Figure  4.17  was  investigated.  The  intratidal  WQM  simulation  was  also 
run  with  a  300  second  time  step  (same  time  step  size  as  the  HM) .  The 
change  of  time  step  did  not  improve  comparisons  with  the  HM.  Thus,  am¬ 
plitude  errors  (numerical  dampening)  associated  with  the  first-order 
(with  time)  Euler  implicit  method  for  WQM  vertical  advection  was  elimi¬ 
nated  as  a  cause  of  model  differences.  However,  the  two  solution 
schemes  can  yield  substantially  different  results  for  areas  where  the 
water  column  is  discretized  with  only  two  or  three  layers.  Amplitude 
and  phase  errors  (related  to  grid  resolution  rather  than  time  step)  of 
advection  differencing  schemes  can  be  quite  large  when  the  physical  wave 
length  being  advectea  is  represented  with  only  a  few  grid  points  (Ander¬ 
son  1984) .  The  physical  wave  in  this  case  is  the  vertical  salinity 
gradient  which  is  advected  up  and  down  by  the  vertical  velocities. 
Therefore,  differences  from  observed  data  and  in  the  two  models  (i.e.  HM 
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and  WQM)  should  not  be  surprising  considering  the  poor  grid  resolution 
in  some  areas  (e.g.  Stations  MP  and  UY) .  A  minimum  of  three  layers 
should  probably  have  been  used  throughout  the  model. 

4.3.3  Benefits  of  Intertidal  Transport 

One  of  the  biggest  advantages  of  using  intertidal,  residual  hydro¬ 
dynamics  for  use  in  a  transport  model  rather  than  using  intratidal  hy¬ 
drodynamics  is  the  reduction  in  computer  hard  disk  space  requirements. 
The  intertidal  hydrodynamics  for  the  1985  simulation,  in  binary  form, 
required  135  megabytes  (MB)  of  disk  space.  Intratidal  hydrodynamics  for 
1985  required  ten  times  as  much  space,  or  1,350  MB  (i.e.  1.35  gigabytes, 
GB) .  Presently,  1.35  GB  of  storage,  even  at  a  supercomputer  site,  is 
considerable.  Applications  with  more  grid  cells  could  result  in 
unacceptable  disk  space  requirements  for  intratidal  averaging  of  long¬ 
term  HM  output.  Therefore,  savings  in  disk  space  are  important. 

Both  the  Eulerian  and  Stokes'  flows  were  written  to  disk,  which  was 
not  efficient  use  of  disk  space.  Both  types  of  flows  were  written  so 
that  they  could  be  compared.  For  production  runs,  the  two  types  of 
flows  could  be  added  together  within  the  processor,  forming  Lagrangian 
residual  flows,  thus,  reducing  disk  storage  space  requirements.  The 
1985  intertidal  simulation  requires  74  MB  with  this  implementation. 

The  CPU  time  savings  on  the  Cray  Y/MP  for  intertidal  WQM  simula¬ 
tions  versus  intratidal  runs  were  not  as  great  as  originally  envisioned. 
For  example,  the  WQM  with  intratidal  (4500  sec)  hydrodynamic  update  per¬ 
iods  required  828  CPU  seconds  versus  488  seconds  with  intertidal  (12.5 
hour)  hydrodynamic  update  periods  for  the  1985  salinity  simulation. 

Thus,  the  intratidal  run  required  about  70  percent  more  CPU  time  than 
the  intertidal  run.  The  average  and  maximum  time  steps  were  3182  and 
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8693  sec,  respectively,  for  the  intertidal  run  and  were  2009  and  4500 
sec  for  the  intratidal  run.  The  jifference  in  the  average  time  step  for 
the  two  runs  accounts  for  most  of  the  difference  in  the  CPU  times.  It 
appears  that  the  Cray  Y/MP  reads  binary  information  very  efficiently 
relative  to  other  machines.  The  savings  in  CPU  time  could  be  much 
greater  with  intertidal  than  with  intratidal  updates  on  other  machines 
with  less  efficient  reading  of  input.  Also,  other  applications  with 
different  grid  resolution  could  result  in  much  greater  differences  in 
intratidal  and  intertidal  allowable  time  step  sizes. 

The  results  of  this  study  demonstrate  that  3D,  Lagrangian  residual 
transport  can  be  accomplished  in  a  practical  manner.  It  is  an  accom¬ 
plishment  to  run  a  3D,  intratidal  HM  for  year-long  periods  generating 
hydrodynamic  input  files  for  UQM  transport  that  are  of  reasonable  size. 
The  procedures  described  herein  provide  an  efficient  and  effective  means 
of  indirectly  coupling  models  with  vastly  different  time  scale  require¬ 
ments  for  tidal  systems. 


CHAPTER  5 


SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 


5 . 1  SUMMARY 

The  main  goal  of  this  study,  which  was  to  develop  a  method  for 
computing  3D  Lagrangian  residual  currents  from  an  intratidal  hydrodynam¬ 
ic  model  for  driving  an  intertidal  transport  model,  was  successfully  ac¬ 
complished.  Although  the  nature  of  tide -induced  residual  currents  had 
been  previously  studied,  there  were  no  known  examples  of  the  use  of  in¬ 
tratidal  HM  information  for  developing  3D  Lagrangian  residual  circula¬ 
tion  to  drive  an  intertidal  transport  model.  The  basic  hypothesis  of 
this  work  (i.e.  tide-induced  residual  currents  are  important  for  Chesa¬ 
peake  Bay)  was  found  to  be  true. 

A  3D  HM  that  uses  boundary- fitted  coordinates  in  planform  was  in¬ 
directly  coupled  to  a  WQM,  which  uses  the  integrated  compartment  method. 
The  coupling  of  the  two  models  was  accomplished  through  the  development 
of  an  interface  processor  implemented  within  the  HM.  The  processor  con¬ 
verts  nondimens ional ,  contravariant  velocities  in  transformed  coordi¬ 
nates  to  dimensional,  physical  flows  for  the  WQM.  Conversion  from  con¬ 
travariant  to  physical  components  retains  flows  normal  to  transverse 
grid  lines  as  required  by  the  WQM. 

The  sum  of  the  Eulerian  residual  velocity  and  the  Stokes'  drift  was 
used  as  a  first-order  approximation  for  the  Lagrangian  residual  cur¬ 
rents.  The  Stokes'  drift  approximates  residual  currents  induced  by  the 
nonlinear  interactions  of  the  tidal  currents  and  represents  the  net 
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drift  experienced  by  a  particle  passing  through  a  spatially  varying 
velocity  field  in  an  oscillating  flow.  A  formulation  for  Stokes'  drift 
was  obtained  which  guarantees  mass  conservation.  The  method  was  imple¬ 
mented  within  the  interface  processor  so  that  intertidal  hydrodynamic 
information  could  be  processed  and  output  as  the  HM  is  running.  This 
information  is  subsequently  used  to  drive  intertidal  WQM  transport. 

The  methods  were  first  tested  against  a  2D  (vertical-longitudinal) 
analytical  result  to  ensure  proper  implementation.  Other  than  some 
adjustment  of  eddy  viscosity  to  account  for  the  effects  of  numerical 
dampening,  the  computed  residual  currents  compared  favorably  with  the 
analytical  result,  thus,  confirming  correct  implementation  of  the  proce¬ 
dures.  The  HM  and  WQM  were  then  applied  to  Chesapeake  Bay  for  the  peri¬ 
od  September  1983.  This  simulation  confirmed  proper  linkage  of  the  two 
models  with  close  agreement  for  intratidal  salinity  transport. 

The  methodology  was  more  fully  evaluated  through  an  application  on 
Chesapeake  Bay  for  the  entire  year  1985.  Salinity  observed,  computed  by 
the  HM,  and  computed  by  the  WQM  was  the  basis  for  making  transport  eval¬ 
uations.  Salinity  computed  with  intertidal  WQM  transport  (i.e.  with 
Lagrangian  residuals)  showed  good  agreement  with  observed  salinity  and 
that  resulting  from  intratidal  transport.  The  mean  absolute  difference 
(i.e.  MAE)  in  salinity  resulting  from  intertidal  versus  intratidal  WQM 
transport  was  0.98  ppt. 

5 . 2  CONCLUSIONS 

Application  of  the  HM  and  interface  processor  to  the  2D  analytical 
test  case  of  Ianniello  (1977)  revealed  that  the  methods  were  correctly 
implemented.  However,  the  results  were  affected  by  the  choice  of  the 
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value  for  the  nondimens ional  vertical  scale  factor,  Zr.  The  conserva¬ 
tive  formulation  for  residual  currents  implemented  in  the  processor 
yielded  results  different  from  the  analytical  result  for  the  Eulerian 
residual  and  Stokes'  drift  in  the  surface  layer.  However,  the  two 
components  added  together  gave  Lagrangian  residuals  quite  similar  to  the 
analytical  result. 

The  21  day  simulation  (i.e.  September  1983)  of  Chesapeake  Bay  was  a 
good  test  case  for  testing  interfacing  of  the  HM  and  WQM,  but  the  period 
was  not  long  enough  to  evaluate  residual  transport  results.  The  year¬ 
long  simulation  of  1985  was  an  excellent  test  case  for  transport  evalua¬ 
tions  . 

Intertidal  salinity  transport  results  without  the  Stokes'  drift 
component  clearly  demonstrated  the  need  to  consider  tide -induced  residu¬ 
al  currents  for  intertidal  transport  modeling.  All  results  indicated 
that  the  basic  effect  of  the  Stokes'  flows  is  to  transport  salinity  up- 
estuary.  Neglecting  Stokes'  flows  in  intertidal  transport  would  tend  to 
under-estimate  salt  water  (or  other  ocean  boundary  solutes)  intrusion 
and  overall  circulation.  The  Stokes'  flows  were  found  to  be  about  the 
same  order  of  magnitude  as  the  Eulerian  flows  in  the  Chesapeake  Bay 
tributaries . 

Although  horizontal  diffusion  was  found  to  be  relatively  unimpor¬ 
tant,  salinity  transport  results  were  sensitive  to  vertical  diffusion. 
However,  computed  intertidal  salinity  transport  was  not  as  sensitive  to 
neglecting  vertical  diffusion  as  to  neglecting  Stokes'  flows.  Simple 
time-averages  of  HM  vertical  diffusivities  were  used  for  intertidal  WQM 
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Intratidal  salinity  transport  was  much  more  sensitive  than  inter¬ 
tidal  transport  to  using  the  first-order,  upwind  differencing  scheme  for 
horizontal  advection.  Considerably  more  salinity  was  numerically  dif¬ 
fused  upstream  in  the  tributaries  for  the  intratidal  results.  Higher- 
order  advection  differencing  schemes  (e.g.  QUICKEST)  are  more  important 
for  intratidal  transport  than  for  intertidal  transport  because  of  the 
greater  velocities  involved  in  intratidal  transport. 

Intertidai  transport  results  were  insensitive  to  variations  in  the 
time  step  size;  this  was  attributed  to  the  use  of  the  QUICKEST  scheme 
for  horizontal  advection.  Although  the  vertical  advection  solution 
scheme  of  the  WQM  was  only  second- order  accurate  with  space  and  first- 
order  accurate  with  time,  higher-order  accuracy  for  vertical  advection 
was  not  nearly  as  important  as  for  horizontal  advection  since  the  verti¬ 
cal  velocities  were  much  smaller. 

The  first-order  estimates  for  Lagrangian  residuals  (i.e.  Eulerian 
residual  plus  Stokes'  drift)  resulted  in  some  loss  in  transport  informa¬ 
tion,  especially  in  the  tributaries.  Differences  in  intertidal  versus 
intratidal  salinity  transport  results  were  partially  attributed  to 
second-order  tidal  phase  effects  (i.e.  Lagrangian  drift)  that  can  not  be 
accounted  for  in  this  first-order  method.  Other  discrepancies  in  in¬ 
tratidal  HM  versus  intertidal  WQM  results,  such  as  in  the  middle  and 
upper  tributaries,  were  attributed  to  differences  in  the  vertical  advec¬ 
tion  differencing  schemes  of  the  HM  and  WQM.  Poor  vertical  grid  resolu¬ 
tion  (i.e.  two  layers)  accentuated  these  differences.  Where  three  or 
more  layers  were  used,  differences  due  to  the  solution  schemes  were  not 
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In  the  Chesapeake  Bay  application,  one  of  the  benefits  of  process¬ 
ing  intratidal  HM  information  for  intertidal  transport  was  the  order  of 
magnitude  reduction  in  disk  file  space  required  to  hold  the  information. 
File  space  requirements  were  reduced  from  approximately  a  gigabyte  to 
about  100  megabytes.  Intratidal  WQM  salinity  simulations  for  Chesapeake 
Bay  required  about  70  percent  more  CPU  time  than  intertidal  WQM  simula¬ 
tions.  Savings  in  CPU  time  was  not  nearly  as  great  as  first  envisioned 
due  to  the  efficient  input-output  capabilities  of  the  Cray  Y/MP.  Great¬ 
er  savings  might  be  realized  on  other  types  of  computers. 

5.3  RECOMMENDATIONS  FOR  FUTURE  STUDIES 

Although  intertidal  salinity  transport  results  showed  relatively 
good  agreement  with  observed  data  and  intratidal  results,  there  was  some 
disparity,  especially  in  the  tributaries  where  greater  tidal  nonlineari¬ 
ties  exist.  In  the  tributaries  of  Chesapeake  Bay,  k  approaches  the 
range  of  0.15  to  0.20.  The  limits  of  this  first-order  approach  imposed 
by  greater  tidal  nonlinearities  (i.e.  larger  k)  need  to  be  further 
studied.  These  tests  should  be  conducted  for  a  simple  test  case,  such 
as  the  2D  channel  of  Chapter  4.  Intratidal  and  intertidal  salinity  (or 
other  tracer)  transport  comparisons  should  be  made  for  a  range  of  k, 
which  can  be  obtained  by  varying  the  tidal  amplitude  at  the  ocean  boun¬ 
dary. 

The  methods  developed  herein  should  be  applied  to  other  tidal  sys¬ 
tems.  Chesapeake  Bay  is  a  partially  mixed  estuary.  Future  applications 
might  include  a  system  that  is  more  stratified  and  one  that  is  more 
fully  mixed.  More  highly  stratified  systems  result  when  the  tidal  cur¬ 
rents  have  less  influence  on  mixing.  Thus,  highly  stratified  estuaries 
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would  be  less  nonlinear  than  Chesapeake  Bay,  and  the  Stokes'  flows 
should  have  less  importance  in  those  systems.  In  a  more  fully  mixed 
estuary,  the  Stokes'  flows  may  not  adequately  describe  the  tide-induced 
residual  transport.  In  those  cases,  intratidal  transport  may  be  neces¬ 
sary. 

Future  work  might  consider  trying  to  develop  a  method  of  processing 
tidally-averaged  hydrodynamics  that  include  second-order,  tidal  phase 
effects  for  application  to  tidal  systems  with  greater  nonlinearity. 
Although  Lagrangian  drift  terms  are  tidal  phase  dependent,  there  may  be 
a  way  to  convert  Lagrangian  drift  into  tidally-averaged  tidal  phase 
dispersion. 

More  work  should  be  done  to  determine  the  effects  of  grid  resol¬ 
ution  on  intertidal  versus  intratidal  numerical  transport  results. 

There  is  a  question  of  whether  the  grid  resolution  affects  the  inter¬ 
tidal,  Lagrangian  residual  transport  differently  from  intratidal  trans¬ 
port.  This  question  can  be  answered  by  comparing  intratidal  versus 
intertidal  transport  for  various  vertical  and  horizontal  grid  densities 
applied  to  simple  geometries. 

Some  three-dimensional  hydrodynamic  models  use  a  grid  transforma¬ 
tion  in  the  vertical  dimension  (i.e.  sigma  stretching).  With  sigma 
stretching,  the  grid  has  the  same  number  of  layers  for  all  planform 
cells,  and  all  layers  expand  and  contract  as  the  water  surface  rises  and 
falls.  Such  a  concept  can  be  considered  a  pseudo- Eulerian- Lagrangian 
grid  scheme  since  the  layers  move  up  and  down  with  the  flow.  It  is 
questionable  whether  the  methods  presented  herein  will  work  with  such  a 
grid.  The  Stokes'  flows  are  based  on  Eulerian  field  velocities.  Be¬ 
cause  of  the  moving  layer  interfaces  in  sigma  stretching,  the  frame  of 
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reference  is  not  fixed.  Therefore,  vertical  velocities  are  not  truly 
Eulerian  field  variables  in  sigma  stretched  coordinates.  A  Lagrangian 
residual  processor  for  a  sigma  grid  HM  should  be  developed  and  tested 
against  the  2D  analytical  result  of  Chapter  4.  The  Lagrangian  residual 
currents  must  be  processed  in  a  manner  that  conserves  mass  to  be  useful 
in  intertidal  transport  computations. 

The  Chesapeake  Bay  model  will  be  used  to  provide  information  for 
nutrient  management.  It  is  recommended  that  hydro - environmental  models, 
such  as  the  Chesapeake  Bay  model,  also  be  used  for  developing  monitoring 
strategies.  Model  studies  usually  follow  monitoring  studies  since  field 
data  is  needed  for  model  calibration.  After  the  models  are  calibrated, 
they  could  be  used  to  develop  more  soundly  based  future  monitoring  stud¬ 
ies.  The  latter  activity  rarely  occurs  now.  Models  can  help  to  better 
understand  the  system,  define  data  gaps  and  needs,  and  determine  impor¬ 
tant  system  features.  The  Chesapeake  Bay  model  and  other  similar  models 
should  be  used  to  guide  future  water  quality  monitoring  efforts. 
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APPENDIX  A 


HYDRODYNAMIC  MODEL  COORDINATE 
TRANSFORMATION  AND  NONDIMENSIONALIZATION 

Basic  concepts  of  the  coordinate  transformations  and  the  nondimen- 
sionalization  for  the  boundary- fitted  hydrodynamic  model,  CH3D,  are 
needed  for  conversion  from  hydrodynamic  model  nondimens ional  quantities 
in  transformed  coordinates  to  dimensional,  physical  quantities  required 
for  the  water  quality  model.  The  details  of  tensor  analysis  and  gen¬ 
eralized  coordinate  transformations  can  be  found  in  Sokolnikoff  (1960) 
and  McConnell  (1957).  Application  of  coordinate  transformations  to  CH3D 
can  be  found  in  Sheng  (1986a  and  1986b)  and  Sheng  and  Hirsh  (1984) . 

Equations  in  transformed  coordinates  can  be  obtained  in  terms  of 
the  contravariant,  covariant,  or  physical  velocity  components  through 
tensor  transformations.  The  contravariant  and  physical  components  are 
locally  orthogonal  to  the  grid  lines,  whereas  the  covariant  components 
generally  are  not.  The  three  components  are  identical  in  a  Cartesian 
coordinate  system.  In  a  general,  non-Cartesian  system,  the  three  compo¬ 
nents  are  expressed  in  terms  of  each  other  as 
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where 

U1  -  contravariant  velocities 
UA  -  covariant  velocties 
U(i)  -  physical  velocities 

There  is  no  summation  on  i  in  Equations  A.l  and  A. 2.  The  metric  tensor, 
g^j ,  is  defined  for  the  two-dimensional  case  of  interest  here  as 
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+  A  j 

(A. 3) 


where  x,y  are  Cartesian  coordinates  and  £ ,r)  are  coordinates  of  the 
transformed  plane.  The  physical  components  of  a  vector  are  the  projec¬ 
tions  of  the  vector  on  the  tangents  to  the  coordinate  curves.  From 
Equation  A.l,  physical  grid  velocities  that  are  locally  orthogonal  to 
the  physical  grid  lines  can  be  obtained  from  the  transformed  contra¬ 
variant  velocities,  which  are  the  dependent  variables  within  CH3D. 

The  determinant  of  the  metric  tensor  (Equation  A. 3)  is 

so  -  ( x«  y*  f  <A-4> 

The  Jacobian  of  the  tranformation,  which  is  gQ^/2j  scales  a  contra¬ 
variant  velocity  to  a  physical  flow  (per  unit  depth)  between  two  grid 
lines 

-  4§o  u  (A-5) 

q^-^V  (A.6) 
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where 

-  physical  flow  per  unit  depth  in  the  £  direction 

-  physical  flow  per  unit  depth  in  the  r?  direction 

U  -  contravariant  velocity  in  the  £  direction 

V  -  contravariant  velocity  in  the  r?  direction 


Equation  A. 5  and  A. 6  are  apparent  from  the  continuity  equation  in  gen¬ 
eral  coordinates.  The  Jacobian  is  also  used  to  obtain  the  surface  area 
within  a  grid  cell 

As  -  ^  d£  dr?  (A. 7) 

where 

d£  -  distance  between  £  grid  lines  (df  “  1) 
dr?  °  distance  between  r?  grid  lines  (dr?  *=  1) 


The  conservation  equations  that  are  solved  within  CH3D  are  in  non- 
dimensional  form.  The  independent  variables  (i.e.  spatial  coordinates 
and  time)  are  made  nondimens ional  as  follows: 


(A. 8) 
(A. 9) 
(A. 10) 
(A. 11) 


where  the  astericks  denote  nondimens ional  quantities,  the  variables  with 
subscript  r  arc  reference  quantities ,  and  f  is  the  Coriolis  parameter  as 
defined  in  Chapter  3.  Nondimens ional  dependent  variables  (i.e.  three 
velocity  components  and  water  surface  displacement)  of  interest  here 


are: 
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u •  _  JL 
ue 

(A. 12) 

v*  ~  X 

Ur 

(A. 13) 

Xr 

-  w  - 

u  z 

r  r 

(A. 14) 

"r  (nfc) 

(A. 15) 

where  Ur  is  a  reference  velocity  and  g  is  the  acceleration  due  to  grav¬ 
ity.  Turbulent  eddy  viscocity  and  diffusivity  are  also  made  nondimen- 
sional  by  dividing  by  reference  values.  Details  of  the  nondiinenslonal 
variables  and  equations  can  be  found  in  Sheng  (1986a  and  1986b)  and 
Johnson  et  al.  (1989). 


APPENDIX  B 


PROCESSOR  PROGRAM  LISTING 

The  processor  consists  of  the  two  subroutines  WQMOUT  and  PROCZ  which 
are  listed  below.  The  subroutines  are  written  in  Fortran  77  and  are 
compiled  separately  and  linked  with  the  compiled  object  file  for  CH3D. 
The  resulting  executable  has  been  run  on  a  VAX  8800  and  Cray  2  and  Cray 
Y/MP  supercomputers  in  batch  mode. 


************************************************************************ 

*  SUBROUTINE  WQMOUT  * 

********* *************************************************************** 

*  * 

*  This  subroutine  takes  CH3D  information  and  processes  it  into  ICM  WQ  * 

*  model  information.  Time- invariant  grid  information  and  time -varying* 

*  transport  information  are  output.  Time-varying  information  consists* 

*  of  3D  Eulerian  and  Stokes'  flows  and  time-averaged  vertical  * 

*  diffusivities .  Intratidal  and  intertidal  averaging  is  done  „ith  * 

*  the  same  subroutine  by  specifying  the  averaging  interval  in  file94.  * 

*  This  subroutine  is  intended  for  use  with  the  Z-grid  (varying  number  * 

*  of  layers)  version  of  CH3D.  This  version  was  written  with  rectan-  * 

*  gular  grid  overlays  in  mind,  but  it  has  not  been  tested  for  overlays* 
************************************************************************ 


*********************************************************************** 


*  GRID  MAPPING  VARIABLES  * 

*  * 

*  NBP  -  max  dimension  on  number  of  boxes  in  surface  layer  * 

*  NFP  -  max  dimension  on  number  of  horizontal  flow  faces  in  * 

*  surface  layer  * 

*  NSR  -  number  of  surface  boxes  * 

*  NB  individual  box  number  in  horizontal  plane  * 

*  IFIRST  -  first  I  index  in  the  hydrodynamic  grid  to  overlay  in  * 

*  box  NB  * 

*  ILAST  -  last  I  index  in  the  hydrodynamic  grid  to  overlay  in  * 

*  box  NB  * 

*  JFIRST  -  first  J  etc  * 

*  JLAST  -  last  J  etc  * 

+  NHQF  -  number  of  horizontal  flow  faces  in  the  surface  layer  * 
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* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


NHQFT 

FN 

QD 

IQ 

JQ 

IL 

JR 

KP 

KZ 

KF,KL 


SB 

K 

F 


-  total  number  of  horizontal  flow  faces 

-  cell  face  number  in  horizontal  plane 

-  flow  direction  code  (X-l ,Y-2 , Z-3) 

-  I  in  box  model  -  box  flow  is  from 

-  J  in  box  model  -  box  flow  is  to 

-  box  model  cell  to  the  left  of  IQ,  when  counting  from 
left  to  right 

-  box  model  cell  to  the  right  of  JQ,  etc 

-  index  of  CH3D  grid  line  perpendicular  to  box  model 
flow  Q(FN) 

-  index  of  CH3D  layer  corresponding  to  box  model  cell 
-*  range  of  CH3D  horizontal  cells  for  spatial  averaging 

of  flows  into  the  box  model.  For  a  1:1  overlay  of 
flows  to  box  model  cells,  KF  -  KL; 
for  QD  -  1 

index 
index 


KP  -  I 
KF,KL  -  J 
for  QD  2 
KP  -  J 
KF,KL  -  I 


index 
index 

counter  designating  surface  box  number 
counter  designating  layer  number 
counter  designating  face  number 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DT 

IM 

JM 

1 1 J  ,K 
KMAX 

KM(I , J) 

XREF 

SREF 

ZREF 

UREF 

AVR 

RB 

U(I,J,K) 

V(I,J,K) 

W(I , J ,K) 

S(I,J) 

ASA(I , J) 
DELTAZ 


HYDRODYNAMIC  MODEL  VARIABLES 
dimensionless  HM  time  step 

number  of  grid  cells  in  the  x(xi)  -  direction  +  1 
number  of  grid  cells  in  the  y(eta)  -  direction  +  1 
CH3D  grid  cell  indices 

maximum  number  of  grid  cells  in  the  z  -  direction 
(also  KMAX  is  the  surface  layer) 

layer  number  for  bottom  layer  of  cell  I , J ;  KM~1  for 
deepest  area 

scale  factor  to  nondimensionalize  x  and  y 
dimensions 

scale  factor  to  nondimensionalize  the  water 
surface  elevation 

scale  factor  to  nondimensionalize  the  z  dimension 
scale  factor  to  nondimensionalize  the  velocities 
scale  factor  for  nondimen  vertical  eddy  diff 
Rossby  number  used  for  nondimensional  scaling 
dimensionless,  contravariant  velocity  on  left  face 
of  cell  I , J ,K  in  planar  (x-dir  vel) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


-  dimensionless,  contravariant  velocity  on  bottom  face* 
of  cell  I , J , K  in  planar  (y-dir  vel)  'n* 

«■  dimensionless,  contravariant  velocity  on  top  face  of* 

* 
* 
* 
* 
* 
* 


cell  I,J,K  in  vertical  plane  (z-dir  vel) 

water  surface  displacement  from  top  of  surface  cell 

I,J 

time-averaged  S 

layer  thickness  for  all  cells  below  the  surface 
layer 
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*  DELTAZM  -  nominal  layer  thickness  for  the  surface  layer  * 

*  DZUU(I,J,K)  -  surface  layer  thickness  located  (averaged)  at  the  * 

*  U  flow  face  * 

*  DZW(I,J,K)  =  surface  layer  thickness  located  (averaged)  at  the  * 

*  V  flow  face  * 

*  GB(I,J,K)  =  dimensionless  vertical  diffusivity  on  top  face  of  * 

*  cell  I , J ,K  * 

*  Gll(I,J,L)=  metric  coefficient  which  scales  the  transformation  * 

*  of  grid  in  the  xi  direction  * 

*  G22(I,J,L)-  metric  coefficient  which  scales  the  transformation  * 

*  of  grid  in  the  eta  direction  * 

*  GD(I,J,L)  -  Jacobian  of  the  metric  tensor  for  grid  transf.  * 

*  * 

*  Note  that  the  L°l,2,3  in  the  metric  coefficient  arrays  specifies* 

*  the  position  on  the  computational  cell  in  which  each  component  * 

*  is  defined.  Specifically:  * 

■*  * 

*  1  -  cell  center  * 

*  2  -  left  cell  face  in  plan  * 

*  3  =  bottom  cell  face  in  plan  * 

*  * 

*  * 

*  PROCESSOR  VARIABLES  * 

*  * 

*  ITWQS  °  time  iteration  when  time-varying  processing  is  ,-.o  * 

*  begin  * 

*  IKNT  **  counter  for  checking  when  the  end  of  the  averaging  * 

*  has  been  reached  * 

*  NAVG  »  number  of  HM  time  steps  to  average  over  * 

*  BL(SB,2)  *=  box  lengths  in  horizontal  directions  for  box  model;  * 

*  Z  direction  will  be  computed  in  WQM  from  CVOLS  and  * 

*  BSAREA  * 


*  CVOL(SB)  =  volume  of  cells  below  the  surface  layer 

*  CVOLS (SB)=  volume  of  cells  in  the  surface  layer 

*  BSAREA(SB)  =■  time  invariant  surface  areas  of  boxes 

*  FAREA(F)  «=■  horizontal  flux  face  area  of  cells 

*  HQ  **  dimensional,  physical  horizontal  flow  for  face  F 

*  ZQ(SB,K)  =  dimensional,  physical  vertical  flow  for  cell 

*  AVGQ(F)  <=  time  average  of  HQ 

*  AVGZQ(SB.K)  =  time  averaged  ZQ 

*  AVDIFZ  "  =  time -averaged,  dimensional  vertical  diffusivity 

*  NCP  =  number  of  corner  points  to  be  read  in  to  zero  Az's 

*  NEX  =  number  of  bottom  corners  to  be  read  in  to  zero  Ax's 

*  NNY  =  number  of  bottom  corners  to  be  read  in  to  zero  Ay's 

*  IC,JC,KC  =  I,J,K  indices  for  Az  corner  points 

*  IE,JE,KE  =  I,J,K  indices  for  Ax  corner  points 

*  IN,JN,KN  -  I,J,K  indices  for  Ay  corner  points 

*  UA(I,J,K)=  time - averaged  U 

*  VA  "  time -averaged  V 

*  WA  "  =  time -averaged  W 

*  UD  "  =  cumulative  displacement  resulting  from  U  velocity 

*  VD  "  =  cumulative  displacement  resulting  from  V  velocity 

*  WD  "  =  cumulative  displacement  resulting  from  W  velocity 

*  UDA  "  =  time -averaged  UD 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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*  VDA  "  -  time-averaged  VD  * 

*  WDA  "  -  time -averaged  WD  * 

*  AX(I,J,K)«*  vector  potential  Ax  for  cell  I,J,K  * 

*  AY(I,J,K)-  vector  potential  Ay  for  cell  I,J,K  * 

*  AZ(I,J,K)-*  vector  potential  Az  for  cell  I,J,K  * 

*  ST(F)  -  horizontal  Stokes  flow  on  face  F  * 

*  STZ(SB,K)-  vertical  Stokes  flow  * 

*  * 


*********************************************************************** 

SUBROUTINE  WQMOUT 

INCLUDE  ' chesv . inc ' 

INCLUDE  'comm3dv. inc' 

COMMON  /AVG/  UA(0:IM,0:JM,0:KMAX) ,  VA(0:IM,0:JM,0:KMAX) , 
WA(0:IM,0:JM,0:KMAX) , 

AX(0:IM,0:JM,0:KMAX) ,  AY(0:IM,0:JM,0:KMAX) , 
AZ(0:IM,0:JM  0:KMAX), 

UD(0:IM,0:JM,0:KMAX) ,  VD(0:IM,0:JM,0:KMAX) , 
WD(0:IM,0:JM,0:KMAX) , 

UDA(0:IM,0:JM,0:KMAX) ,  VDA(0:IM,0:JM,0:KMAX) , 
WDA(0:IM,0:JM,0:KMAX) ,  ASA(0:IM,0:JM) , 

NAVG,  NCP, 

NEX,  NNY, 

IC(500),  JC(500) ,  KC(500) 

IE(500),  JE(500) ,  .  KE(500'> 

IN(500),  JN(500) ,  KN(50f 


PARAMETER  (NBP-800 .NFP-1500) 

DIMENSION  AVGQ(NFP*KMAX) ,  AVGZQ(NBP,KMAX) ,  AVDIFZ(NBP.KMAX) 

INTEGER  QD,  F,  FN,  SB,  TVD 
DIMENSION  BL(NBP,2),  CVOL(NBP) , 

CVOLS(NBP) ,  BSAREA(NBP) ,  ZQ(NBP,KMAX) , 

ST(NFP*KMAX)  ,  STZ(NBP.KMAX) 

DIMENSION  IQ(NFP*KMAX) , JQ(NFP*KMAX) ,  IL(NFP*KMAX) , 
JR(NFP*KMAX) 

COMMON  /B0X01/  IFIRST(NBP),  ILAST(NBP) ,  JFIRST(NBP) , 

JLAST(NBP) ,  NSB 

COMMON  /B0X02/  NHQF,  QD(NFP*KMAX) ,  KP(NFP*KMAX) , 

KF (NFP*KMAX) ,  KL(NFP*KMAX) ,  KZ(NFP*KMAX) 

COMMON  /B0X03/  SRG11(0:IM,0:JM,3) ,  SRG22(0:IM,0:JM,3) 

COMMON  /BOX04/  AVGN,  AVGQ, 

FAREA(NFP*KMAX) ,  AVGZQ , 

AVDIFZ 


DATA  TVD  /96/ 
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*********************************************************************** 

*  TASK  1:  INPUT  SECTION  * 

*********************************************************************** 
C 

C***  READ  IN  CORNER  POINTS  FOR  A  TERMS 
C 

READ(85 ,*)  NCP 
DO  30  N=1 ,NCP 

READ(85 ,*)  NP, IC(N) , JC(N) ,KC(N) 

30  CONTINUE 

READ(85 ,*)  NEX 
DO  32  N-l.NEX 

READ(85 ,*)  NP,IE(N) ,JE(N) ,KE(N) 

32  CONTINUE 

READ(85 ,*)  NNY 
DO  34  N=1 , NNY 

READ(85 ,*)  NP,IN(N) ,JN(N) ,KN(N) 

34  CONTINUE 

READ (94,*)  NSB,  NAVG ,  ITWQS 
DO  40  SB=1 ,NSB 

READ(94 ,*)  NB, IFIRST(SB) , ILAST(SB) , JFIRST(SB) ,JLAST(SB) 

40  CONTINUE 
READ(95 ,45) 

45  FORMAT(5(/) ) 

READ (95,*)  NHQF , NHQFT 
READ(95 ,46) 

46  FORMAT ( 8 0A1) 

DO  50  F“1 , NHQFT 

READ(95 ,*)  FN,QD(F) , IL(F) , IQ(F) , JQ(F) , JR(F) , 

KP(F) ,KF(F) ,KL(F) ,KZ(F) 

50  CONTINUE 

*********************************************************************** 

*  TASK  2:  INITIALIZATION  SECTION  * 

*********************************************************************** 

AVGN  -  FLOAT (NAVG) 

C 

C***  ZERO  COMPUTATIONAL  VARIABLES 
C 

IKNT  »  0 
DO  9000  1=1, IM 
DO  9000  J=1 , JM 

IF(KM(I , J)  .NE.  0)  THEN 
AX(I, J , KM(I , J ) - 1)=0 . 0 
AY(I,J,KM(I,J)-1)=0.0 
END  IF 

DO  8000  K-l.KMAX 
UA(I , J ,K)=0 . 0 
VA(I,J,K)“0.0 
WA(I ,  J  ,K)**0 .0 
AX(I , J ,K)=0 . 0 
AY ( I , J , K) =0 . 0 
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AZ(I , J ,K)-0 . 0 
UD(I,J ,K)-0.0 
VD(I , J ,K)-0 . 0 
WD(I,J,K)-0.0 
UDA(I , J ,K)-0 . 0 
VDA(I,J ,K)-0.0 
WDA(I , J ,K)-0 . 0 
8000  CONTINUE 
9000  CONTINUE 

DO  10010  F-l , NHQFT 
AVGQ(F)  -  0.0 
10010  CONTINUE 

DO  10030  SB-1, NSB 

DO  10020  F-KM(IFIRST(SB) , JFIRST(SB) ) ,KMAX 
AVGZQ(SB.F)  -  0.0 
AVDIFZ(SB.F)  -  0.0 
10020  CONTINUE 
10030  CONTINUE 

*********************************************************************** 

*  TASK  3:  TIME -INVARIANT  CALCULATIONS  * 

*********************************************************************** 

C 

C***  COMPUTE  X-Y  SCALING  FACTORS 
C 

DO  10060  I-l.IM-l 
DO  10050  J-l.JM-l 
DO  10040  K-1,3 

SRG11(I , J ,K)  -  SQRT(G11(I , J ,K) ) 

SRG22(I , J ,K)  -  SQRT(G22(I , J ,K) ) 

10040  CONTINUE 
10050  CONTINUE 
10060  CONTINUE 

C 

C***  COMPUTE  CELL  SURFACE  AREAS  AND  CELL  LENGTHS 
C 

DO  10130  SB-1, NSB 
I-IFIRST(SB) 

J-JFIRST(SB) 

BSAREA(SB)  =  GD ( I , J , 1 ) *XREF*XREF 
C 

C***  COMPUTE  X- DIRECTION  BOX  LENGTHS 
C 

BL(SB,1)  -  XREF*SRG11(I,J,1) 

C 

C***  COMPUTE  Y- DIRECTION  BOX  LENGTHS 
C 

BL(SB ,2)  -  XREF*SRG22 ( I , J , 1) 

10130  CONTINUE 
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C 

C***  WRITE  INITIAL  DATA  TO  DISK 
C 


WRITE(TVD) 
WRITE (TVD) 
WRITE(TVD) 
WRITE(TVD) 


DELTAZ*0 . 01 

(BSAREA(SB)*1 . 0E-04 , SB-1 ,NSB) 
(BL(SB,1)*. 01, SB-1, NSB) 
(BL(SB, 2)* . 01 , SB-1, NSB) 


C***  FORMATTED  WRITES 


C 

C 

C 

C 


WRITE(88 , 3000) 
WRITE(88 , 3050) 
WRITE(88 , 3000) 
WRITE(88 , 3000) 


DELTAZ*0 . 01 

(BSAREA(SB)*1.0E- 04, SB-1, NSB) 
(BL(SB,1)*. 01, SB-1, NSB) 

(BL(SB , 2)*.01 , SB-1 , NSB) 


300  FORMAT (IX, All/ (IX, 10 (1PE13 . 6)) ) 

305  FORMAT ( IX , All/ ( IX , 8 ( 1PE15 . 8 ) ) ) 

3000  FORMAT(10(1PE13 . 6)  ) 

3050  FORMAT(8(lPE15 . 8) ) 


RETURN 


*********************************************************************** 


*  ENTRY  WQMCVOL  * 
C***  COMPUTE  CVOL,  CVOLS,  AND  FAREA  FIRST  TIME  * 
*********************************************************************** 


ENTRY  WQMCVOL 


C 

C***  CALCULATE  HORIZONTAL  FLUX  FACIAL  AREAS 
C 

DO  10133  F-l.NHQFT 
KSS  =  KF(F) 

IF  (QD(F).EQ.l)  THEN 

FAREA (F)  -  DZUU(KP(F) ,KSS ,KZ(F) )*ZREF*SRG22(KP(F) ,KSS , 2)*XREF 
ELSE 

FAREA (F)  -  DZW(KSS,KP(F),KZ(F))*ZREF*SRG11(KSS,KP(F),3)*XREF 
END  IF 
10133  CONTINUE 

C 

C***  COMPUTE  CELL  VOLUMES 

n 


DO  11130  SB-1, NSB 
I-IFIRST(SB) 

J-J FIRST (SB) 

CVOLS (SB)  -  (DELTAZM  +  S(I,J)*SREF) 

*  GD ( I , J , 1 ) *XREF*XREF 
CVOL (SB)  -  DELTAZ*GD ( I , J , 1 ) *XREF*XREF 
11130  CONTINUE 
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WRITE (TVD)  (FAREA(F)*1 . 0E-04 , F-l , NHQFT) 

WRITE(TVD)  (CV0L(SB)*1 . 0E-06 , SB-1 ,NSB) 

WRITE (TVD)  (CVOLS (SB)*1 . 0E-06 , SB-1 ,NSB) 

C***  FORMATTED  WRITES 

C  WRITE(88 , 3050)  (FAREA(F)*1.0E-4,F-1,NHQFT) 

C  WRITE(88,3050)  (CVOL(SB)*l . 0E-06 , SB-1 ,NSB) 

C  WRITE(88 , 3050)  (CVOLS(SB)*1.0E-06 ,SB-1,NSB) 

RETURN 

*********************************************************************** 
*  ENTRY  WQTVD  -  TIME -VARYING  CALCULATIONS  * 

*********************************************************************** 

ENTRY  WQTVD 

IKNT  -  IKNT  +  1 

FACT  -  UREF  *  ZREF  *  XREF  *  RB 


C 

C***  UPDATE  AVERAGES  AND  DISPLACEMENT  INTEGRAL 
C 


10134 

10135 


DO  10135  I 
DO  10135  J 
ASA(I.J) 
DO  10134  K 
UA(I , J ,K) 
VA(I,J,K) 
WA(I,J,K) 
UD(I , J ,K) 
VD(I , J ,K) 
WD(I , J ,K) 
UDA(I , J ,K) 
VDA(I , J ,K) 
WDA(I , J ,K) 
CONTINUE 
CONTINUE 


-  1,  IM 

-  1 ,  JM 

-  ASA(I.J)  +  S (I , J)/AVGN 

-  KM(I , J) ,KMAX 

-  UA(I , J ,K)  +  U(I , J ,K)/AVGN 
+ 

+ 

+ 


VA(I,J,K) 

WA(I,J,K) 

UD(I,J,K) 


V(I , J ,K)/AVGN 
W(I , J ,K)/AVGN 
U(I , J ,K)*DT/2 . 0 

-  VD(I , J ,K)  +  V(I,J,K)*DT/2.0 

-  WD(I , J ,K)  +  W(I , J ,K)*DT/2 .0 
-UDA(I , J ,K)  +UD(I , J ,K)/AVGN 
-VDA(I , J ,K)  +VD ( I , J , K) /AVGN 
-WDA(I , J ,K)  +WD ( I , J , K) /AVGN 


CALL  PROCZ 


C 

C***  COMPLETE  UPDATE  OF  DISPLACEMENT  INTEGRATION 
C 


DO  10138  I  -  1 , IM 
DO  10138  J  -  1 , JM 
DO  10138  K  -  KM(I , J) ,KMAX 


UD(I , J ,K)  -  UD(I , J ,K)  +  U(I,J,K)*DT/2.0 
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VD(I , J ,K)  -  VD(I , J ,K)  +  V(I , J ,K)*DT/2 .0 
WD(I,J,K)  -  WD(I,J,K)  +  W(I,J,K)*DT/2.0 


10138  CONTINUE 
C 

C***  COMPLETE  FINAL  UPDATE  AND  SUBTRACT  TIDAL  AVERAGE  ESTIMATE 
C 


IF(IKNT  .EQ.  NAVG)  CALL  PROC2 


C 

C***  COMPUTE  HORIZONTAL  FLOWS 
C 


DO  10230  F  -  l.NHQFT 


C 

C***  X-DIRECTION 
C 

IF  (QD(F).EQ.l)  THEN 
KSS  -  KF(F) 

IF(F  .GT.  NHQF)  THEN 
ZSCL  -  DELTAZ 
ELSE 

ZSCL  -  DZUU(KP(F) ,KSS ,KZ(F) )*ZREF 
END  IF 

HQ  -  U(KP(F) ,KSS ,KZ(F) )*UREF*GD(KP(F) ,KSS , 2)* 
ZSCL*XREF 

AVGQ(F)  -  AVGQ(F)  +  HQ/AVGN 
C 

C***  X  STOKES  DRIFT 
C 


IF(IKNT  .EQ.  NAVG)  THEN 
IF(F  .GT.  NHQF)  THEN 
ZSCLN  =  DELTAZ/ZREF 
ELSE 

ZSCLN  =1.0 
END  IF 


om  /  r?\ 


((  AZ(KP(F)  ,KS.°  .  '  KZ(F))-AZ(KP(F),KSS, 
*  ZSCLN  -  *Y(KP(F),KSS,KZ(F))  + 
AY(KP(F) ,KSS ,KZ(F) -1) )  *  FACT 


W7  fP\  N  N 
\l!  /  /  / 


END  IF 


C 

C***  Y- DIRECTION 
C 


256 


ELSE 

KSS-KF(F) 

IF(F  .GT.  NHQF)  THEN 
ZSCL  -  DELTAZ 
ELSE 

ZSCL  -  DZW(KSS ,KP(F)  ,KZ(F))*ZREF 
END  IF 

HQ  -  V(KSS,KP(F) ,KZ(F))*UREF*GD(KSS,KP(F) ,3)* 
ZSCL*XREF 

AVGQ(F)  -  AVGQ(F)  +  HQ/AVGN 


C 

C***  Y  STOKES  DRIFT 
C 


IF(IKNT  .EQ.  NAVG)  THEN 
IF(F  .GT.  NHQF)  THEN 
ZSCLN  -  DELTAZ/ZREF 
ELSE 

ZSCLN  -  1.0 
END  IF 

ST(F)  -  (  AX(KSS,KP(F) ,KZ(F))  -  AX(KSS ,KP(F) ,KZ(F) -1) 

-  (  AZ(KSS+1 ,KP(F) ,KZ(F))  -AZ(KSS ,KP(F) ,KZ(F) ) ) 
*  ZSCLN)  *  FACT 


END  IF 


END  IF 

10230  CONTINUE 
C 

C***  COMPUTE  VERTICAL  FLOWS  AND  VERTICAL  D1FFUSIVITIES 
C 


DO  10270  SB=1,NSB 
I-IFIRST(SB) 

J-JFIRST(SB) 

DO  10260  F=KM(I , J) ,KMAX 

ZQ(SB,F)  -  W(I , J , F)*GD(I , J , 1) 

*XREF*ZREF*UREF 

AVGZQ(SB.F)  -  AVGZQ(SB.F)  +  ZQ(SB,F)/AVGN 
AVDIFZ(SB.F)  -  AVDIFZ(SB.F)  +  GB ( I , J , F) /AVGN*AVR 


C 

C***  Z  STOKES  DRIFT 
C 
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IF(IKNT  .EQ.  NAVG)  THEN 

STZ(SB,F)  -  (  AY(I+1,J,F)  -  AY(I,J,F)  -  AX(I,J+1,F) 
+  AX(I,J,F)  )  *  FACT 

END  IF 

10260  CONTINUE 
10270  CONTINUE 

**************************************************************** 
C***  OUTPUT  BOX  MODEL  INFO  AT  EACH  AVERAGING  INTERVAL  * 

C  * 

**************************************************************** 

IF  (IKNT  .EQ.  NAVG)  THEN 

C 

C***  COMPUTE  SURFACE  VOLUME 
C 

DO  10330  SB-1, NSB 
I-IFIRST(SB) 

J-JFIRST(SB) 

CVOLS(SB)  -  (DELTAZM+S(I , J)*SREF) 

*GD ( I , J , 1 ) *XREF*XREF 

10330  CONTINUE 
C 

C***  COMPUTE  SURFACE  FACE  AREAS 
C 

DO  10340  F-l.NHQF 
KSS  -  KF(F) 

IF  (QD(F).EQ.l)  THEN 

FAREA(F)  -  SRG22(KP(F) ,KSS , 2)*XREF 

*  DZUU(KP(F) ,KSS ,KZ(F) )*ZREF 

ELSE 

FAREA(F)  -  SRG11(KSS ,KP(F) ,3)*XREF 

*  DZW(KSS  ,KP(F)  ,KZ(F) )*ZREF 

END  IF 
10340  CONTINUE 
C 

C***  WRITE  TIME  VARYING  DATA  TO  WATER  QUALITY  MODEL 
C 


WRITE(TVD) 
WRITE (TVD) 
WRITE (TVD) 
WRITE (TVD) 
WRITE (TVD) 


WRITE (TVD) 


TIME 

(FAREA(F)*1 . 0E-04 , F-l ,NHQF) 

(CVOLS(SB)*l. 0E-06 , SB-1 , NSB) 

(AVGQ(F)*1 . 0E-06 , F-l ,NHQFT) 

( (AVDIFZ(SB , F)*l . 0E-04, 

F-KM(IFIRST(SB) , JFIRST(SB) ) .KMAX-l) , 
SB-1, NSB) 

( (AVGZQ(SB, F)*l . 0E-06 , 

F-KM(IFIRST(SB) .JFIRST(SB) ), KMAX-l) , 
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SB-1, NSB) 

WRITE (TVD)  ( (STZ(SB , F)*l .0E-06 , 

F-KM(IFIRST(SB) .JFIRST(SB)) .KMAX-l) , 
SB-l.NSB) 

WRITE(TVD)  (ST(F)*1.0E-06 , F-l ,NHQFT) 

C  FORMATTED  WRITES 

C  WRITE (8 8, 3000) 

C  WRITE(88 , 3050) 

C  WRITE(88 , 3050) 

C  WRITE (8 8, 3000) 

C  WRITE(88 , 3000) 

C 
C 

C  WRITE(88 , 3000) 

C 
C 

C  WRITE(88 , 3000) 

C 
C 

C  WRITE(88 , 3000) 

C 

C***  RESET  TIME  AVERAGES  TO  ZERO 
C 

IKNT  -  0 

DO  10350  F-l.NHQFT 
AVGQ(F)  -  0.0 
10350  CONTINUE 

DO  10370  SB-l.NSB 
I-IFIRST(SB) 

J-JFIRST(SB) 

DO  10360  F-KM( I, J), KMAX-l 
AVGZQ(SB.F)  -  0.0 
AVDIFZ(SB.F)  -  0.0 
10360  CONTINUE 

10370  CONTINUE 

DO  103SQ  I  -  1 , IM 
DO  10390  J  -  1 , JM 
ASA(I,J)  -  0.0 
DO  l0iS0  K  -  KM(I , J) ,KMAX 

*»ta  /  x  t  va  —  n  n 

Wii\  x  ,  u  ,  i\/  —  v  i  v 

VA(I , J ,K)  -  0.0 
WA(I,J,K)  -  C.O 
AX(I , J ,K)  -  0.0 
AY( I , J , K)  »  0.0 
AZ ( I , J , K)  -  0.0 
UD(I , J ,K)  =  0.0 
VD(I , J ,K)  -  0.0 


TIME 

( FAREA ( F) *1 . OE - 4 , F-l , NHQF) 

(CVOLS (SB) *1 . 0E-06 , SB-1 , NSB) 

(AVGQ(F)*1 . 0E-06 ,F-1,NHQFT) 
((AVDIFZ(SB,F)*1.0E-04, 

F-KM(IFIRST(SB) .JFIRST(SB) ), KMAX-l) , 
SB-1, NSB) 

( (AVGZQ(SB , F)*l. 0E-06 , 

F-KM(IFIRST(SB) , JFIRST(SB) ) , KMAX-l) , 
SB-1, NSB) 

((STZ(SB,F)*1.0E-06, 

F-KM(IFIRST(SB) .JFIRST(SB) ), KMAX-l) , 
SB-1, NSB) 

(ST(F)*1.0E-06,F-1,NHQFT) 
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WD(I,J,K)  -  0.0 
UDA(I , J ,K)  -  0.0 
VDA(I , J ,K)  -  0.0 
WDA(I,J,K)  -  0.0 
10380  CONTINUE 
10390  CONTINUE 

END  IF 

RETURN 

END 

************************************************************ 

*  SUBROUTINE  PROCZ  * 

************************************************************ 


SUBROUTINE  PROCZ 

INCLUDE  'chesv.inc' 
INCLUDE  'comm3dv. inc' 


COMMON  /AVG/  UA(0:IM,0:JM,0:KMAX) , 

WA(0: IM,0: JM,0:KMAX) , 
AX(0:IM,0:JM,0:KMAX) , 

AZ (0 : IM , 0 : JM , 0 : KMAX) , 
UD(0:IM,0:JM,0:KMAX) , 

WD ( 0 : IM , 0 : JM , 0 : KMAX) , 
UDA(0:IM,0:JM,0:KMAX) , 
WDA(0:IM,0:JM,0:KMAX) , 

NAVG,  NCP, 

NEX,  NNY, 

IC(500) ,  JC(500) , 

IE(500),  JE(500) , 

IN(500) ,  JN(500) , 


VA(0:IM,0:JM,0:KMAX) , 

AY ( 0 : IM , 0 : JM , 0 : KMAX) , 

VD(0: IM,0: JM,0:KMAX) , 

VDA(0:IM,0: JM,0:KMAX) , 
ASA(0:IM,0:JM) , 


KC(500) , 
KE(500) , 
KN(500) 


MULT  =  4  *  NAVG 


C 

C***  X- DIRECTION  ***  COMPUTE  AY  AND  AZ  COMPONENTS 
C 

DO  150  J  -  1,  J CELLS 
DO  150  1=1,  ICELLS 

IF(HS(I , J)  .EQ.  0.0  .OR.  HS(I-l.J)  .EQ.  0.0)  GO  TO  150 
C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(J  .EQ.  4  .AND.  I  .LE.  5)  GO  TO  150 

C  ****** 


n 

C***  LEVEL  K  ■=  KM  THRU  KMAX  -  1 
C 

DO  120  K  -  KM(I , J) ,  KMAX-1 

UDZA  -  UD(I , J ,K)  +  UD(I,J,K+1) 


o  o  o 
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WXA  -  W(I,J,K)  +  W(I-3  ,J,K) 

C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(I  .EQ.  6  .AND.  J  .EQ.  4)  WXA  -  2.0*W(I,J,K) 

C  ******** 

AY(I,J,K)  -  AY(I,J,K)  +  UDZA*WXA/MULT 

UYA  -  U(I , J ,K)  +  UU.J-l.K) 

VDXA  -  VD(I,J,K)  +  VD(I-1,J,X) 

AZ(I ,  J ,K)  ~  AZ<I,J,K)  +  UYA*VDXA/MULT 

120  CONTINUE 


***  LEVEL  K  -  KMAX,  AY  IS  ZERO 


K  -  KMAX 

UYA-  U(I,J,K)  +  U(I,J-1,K) 

VDXA  -  VD(I,J,K)  +  VD(I-1,J,K) 

AZ(I , J ,K)  -  AZ ( I , J , K)  +  UYA*V DXA/MULT 

150  CONTINUE 

C 

C***  Y- DIRECTION  ***  COMPUTE  AX  COMPONENT 
C 

DO  250  I  -  1,  ICELLS 
DO  250  J  -  1,  J CELLS 

IF(HS(X,J)  .EQ.  0.0  .OR.  HS(I,J-1)  .EQ.  0.0)  GO  TO  250 
C 

C***  LEVEL  K  -  KM  THRU  KMAX  -  1 
C 

DO  220  K  -  KM(I , J) ,  KMAX  -  1 

VZA  -  V(I,J,K)  -rV(I,J, K+l ) 

WDYA  -  WD(I , J ,K)  +  WD(I , J-l ,K) 

C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(J.  ,LE.  5  .AND.  J  .EQ.  5)  WDYA  -  2.0  *  WD(I,J,K) 

C  ******** 

AX(I , J ,K)  -  AX(I,J,K)  +  VZA*WDYA/MULT 
220  CONTINUE 

C 

C***  K  -  KMAX,  AX  IS  ZERO 
C 

250  CONTINUE 


RETURN 
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*************************  ************************************** 

*  ENTRY  PR0C2  * 

*************************************************************** 

ENTRY  PR0C2 

MULT  -  4  *  NAVG 

C 

C***  ACCUMULATE  LAST  TIME  STEP,  SUBTRACT  MEANS  AND  SCALE  BY  THE 
C***  DEPTH  AND  JACOBIAN 
C 

C 

C***  X- DIRECTION  ***  COMPUTE  AY  AND  AZ  COMPONENTS 
C 

DO  J50  J  »  1,  J CELLS 
DO  350  I  -  1,  ICELLS 

IF(HS(I , J)  .EQ.  0.0  .OR.  HS(I-1,J)  .EQ.  0.0)  GO  TO  350 
C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(J  .EQ.  4  .AND.  I  .LE.  5)  GO  TO  350 

C  ****** 

C 

C***  LEVEL  K  ~  KM  THRU  KMAX  -  1 
C 

GAZ  -  (GD(I, J,2)  +  GD(I,J-1,2) 

+  GD ( I , J , 3 )  +  GD(I-l,J,3))/4. 

DO  320  K  =  KM(I , J) ,  KMAX-1 

UOZA  -  UD(I , J , K)  +  UD ( I , J , K+l ) 

WXA  -  W(I,J,K)  +  W(I- 1 , J ,K) 

C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(I  .EQ.  6  .AND.  J  .EQ.  4)  WXA  =  2.0*W(I,J,K) 

C  ********* 

AY(I , J ,K)  -  AY(I,J,K)  +  UDZA*WXA/MULT 

UYA  -  U(I,J,K)  +  U(I,J-1,K) 

VDXA  =  VD(I,J,K)  +  VD(I-1 , J ,K) 

AZ(I,J,K)  -  AZ(I,J,K)  +  UYA*VDXA/MULT 

C 

C***  SPATIAL  AVERAGE  MEAN  QUANTITIES 
C 

WAX  -  WA(I,J,K)  +  WA(I-1,J,K) 

C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(I  .EQ.  6  .AND.  J  .EQ.  4)  WAX  -  2 . 0*WA(T. ,  J ,K) 

C  ******* 

UDAZ  -  UDA(I,J,K)  +  UDA(I,J,K+1) 
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AY(I , J ,K)  -  (AY(I , J ,K)  -  WAX*UDAZ/4.0)*GD(I,J,2) 

UAY  -  UA(I,J,K)  +  UA(I , J-l ,K) 

VDAX  -  VDA(I,J,K)  +  VDA(I-1,J,K) 

AZ(I , J ,K)  -  (AZ(I , J ,K)  -  UAY*VDAX/4.0)*GAZ 

320  CONTINUE 

C 

C***  LEVEL  K  -  KMAX,  AY  IS  ZERO 
C 

K  -  KMAX 

HAZ  -  (DELTAZM  +  0. 25*(ASA(I , J)  +  ASA(I.J-l) 

+  ASA(I-l.J-l)  +  ASA(I-1 , J) )*SREF)/ZREF 
UAY  -  UA(I ,  J  ,K)  +  UA('£ ,  J-l  ,K) 

VDAX  -  VDA(I , J ,K)  +  VDA(I-1,J,K) 

AZ(I,J,K)  =  (AZ(I , J , K)  -  UAY*VDAX/ 4.0) *GAZ*HAZ 

350  CONTINUE 

C 

C***  Y- DIRECTION  ***  COMPUTE  AX 
C 

DO  450  I-l,  ICELLS 
DO  450  J-l,  J  CELLS 

IF(HS(I,J)  .EQ.  0.0  .OR.  HS(I,J-1)  .EQ.  0.0)  GO  TO  450 
C 

C***  LEVEL  K  -  KM(I , J)  THRU  KMAX  -  1 
C 

DO  420  K  -  KM(I,J),  KMAX  -  1 

VZA  -  V(I,J,K)  +  V ( I , J , K+l ) 

WDYA  -  WD(I,J,K)  +  WD(I , J-l ,K) 

C  SPECIFIC  TC  CHESAPEAKE  BAY 

IF(I  .LE.  5  .AND.  J  .EQ.  5)  WDYA  -  2.0  *  WD(I,J,K) 

C  ****** 

AX(I , J ,K)  -  AX(I , J ,K)  +  VZA*WDYA/MULT 

VAZ  -  VA(I,J,K)  +  VA(I,J,K+1) 

WAY  »  WDA(I , J  ,K)  +  WDA(I,J-1,K) 

C  SPECIFIC  TO  CHESAPEAKE  BAY 

IF(I  .LE.  5  .AND.  J  .EQ.  5)  WDAY  =  2.0  *  WDA(I,J,K) 

C  ********* 

AX(I,J,K)  -  (AX(I , J ,K)  -  VAZ*WDAY/4 . 0 ) *GD ( I , J , 3 ) 


o  o  o 
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420  CONTINUE 
450  CONTINUE 
C 

C***  ZERO  RIVER  AND  BOUNDARY  INFLOW  AX,  AY,  AZ 
C 

DO  500  K  -  1  ,  KMAX 
AX(2 , 35 ,K)  -  0.0 

AX(6 , 28 ,K)  »  0.0 

AX(8,28,K)  -  0.0 

AX(17,32,K)  -  0.0 
AX(22 , 35 ,K)  -  0.0 
AX(62 , 15 ,K)  -  0.0 
AX(48 , 17 ,K)  -  0.0 
AX(49 , 17 ,K)  -  0.0 
AX(50 , 17 ,K)  “  0.0 
AX(37 , 1 ,K)  »  0.0 

AY(41 , 15 ,K)  -  0.0 
AZ(1 , 5 ,K)  “  0.0 

AZ(1 , 37 ,K)  -  0.0 

AZ(1 , 38 ,K)  -  0.0 

AZ(2 , 35 ,K)  -  0.0 

AZ(3 , 35 ,K)  -  0.0 

AZ(6,4,K)  »  0.0 

AZ(6,28,K)  -  0.0 

AZ(7 , 28 ,K)  -  0.0 

AZ(8 , 28 ,K)  -  0.0 

AZ(9 , 28 ,K)  =  0.0 

AZ(17 , 32 ,K)  °  0.0 
AZ(18 , 32 ,K)  -  0.0 
AZ(22 , 35 ,K)  =  0.0 
AZ(23 , 35 ,K)  -  0.0 
AZ(48 , 1/ , K)  =  0.0 
AZ(51 , 17 ,K)  -  0.0 
AZ(62 , 15 ,K)  -  0.0 
AZ(63 , 15 ,K)  -  0.0 
AZ(37,1,K)  -  0.0 

AZ(38,1,K)  -  0.0 

AZ(41 , 15 ,K)  -  0.0 
AZ(41 , 16 ,K)  -  0.0 

500  CONTINUE 


AZ' 


n 

O 


DO  400  N  -  l.NCP 

AZ(IC(N) ,JC(N) ,KC(N))  »  0.0 
400  CONTINUE 


POINTS 
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C 

C***  ZERO  EAST-WEST  AX'S  ON  BOTTOM  CORNERS 
C 

DO  600  N  -l.NEX 

AX(IE(N) ,Jfi(N),KE(N))-0.0 
600  CONTINUE 

C 

C***  ZERO  NORTH -SOUTH  AY'S  ON  BOTTOM  CORNERS 
C 

DO  700  N  -l.NNY 

AY(IN(N) , JN(N) ,KN(N) )-0 . 0 
700  CONTINUE 


RETURN 

END 


